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INTRODUCTION 


This  project  is  concerned  with  the  numberical  simulation  of  the  time- 
dependent  behavior  of  spacecraft  sheaths  and  charging  effects.  The  objec¬ 
tives  are  to  develop  a  time-dependent  simulation  model  for  plasma-spacecraft 
interactions,  and  to  apply  this  model  to  a  geometrical  representation  of  the 
SCATHA  spacecraft.  An  appropriate  representation  is  that  of  a  short  right- 
circular  cylinder  (or  "pillbox"),  of  comparable  length  and  diameter  (see 
Fig.  1),  with  azimuthal  syimetry  so  that  the  potential  and  charge  distribu¬ 
tions  can  be  defined  on  a  grid  In  r-z  space. 

A  computer  program  has  been  developed  for  the  study  of  the  time- 
dependent  sheath.  A  grid  Is  used  to  define  the  spatial  distributions  of 
potential  and  charge  density  In  the  space  around  the  spacecraft.  Figure  2 
illustrates  the  nature  of  the  r-z  grid  representation  used.  The  geometry 
is  axially-symmetric,  with  the  axis  shown  as  the  vertical  dotted  boundary 
line  on  the  left,  labelled  "West".  The  boundary  condition  representing  the 
condition  on  the  potential  at  Infinity  is  applied  to  the  other  boundary 
lines  of  the  grid.  The  Inner  boundary  represents  the  satellite  surface,  on 
the  grid  points  of  which  the  surface  potentials  are  defined.  Associated 
with  each  spatial  grid  point  is  a  volume  of  revolution  in  the  shape  of  a 
torus  of  rectangular  cross-section  (shown  as  shaded  boxes  surrounding  some 
of  the  grid  points).  Only  24  grid  points  are  shown  in  Fig.  2,  for  the  pur¬ 
pose  of  clarity.  In  an  actual  problem  many  more  points  are  used. 

An  important  feature  of  this  grid  representation  is  that  the  zoning  is 
nonuniform.  This  allows  for  fine  zoning  in  the  regions  of  interest,  e.g., 
where  there  are  large  gradients,  and  coarse  zoning  elsewhere,  and  has  the 
advantage  of  optimum  use  of  a  given  number  of  grid  points,  and  therefore 
computer  efficiency  in  large  problems. 

The  plasma  electrons  and  ions  are  simulated  by  a  number  of  discrete 
"computer  particles"  injected  through  the  outer  grid  boundaries.  These  par¬ 
ticle  trajectories  are  followed  step  by  step  through  the  grid.  Similarly, 
emitted  electrons  and/or  Ions  are  injected  from  the  inner  (spacecraft)  sur- 
fact.  Each  simulation  particle  represents  a  definite  number  of  real  parti¬ 
cles,  retaining  the  charge-to-mass  ratio  of  the  real  particles.  The  density 
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of  computer  particles  associated  with  any  grid  point  is  determined  at  any 
instant  of  time  by  the  number  of  these  particles  in  the  box  surrounding  that 
grid  point  at  that  time.  The  potential  distribution  is  found  for  the  par¬ 
ticular  charge  distribution  and  boundary  conditions  at  any  given  time  by 
solving  the  discretized  analog  of  the  Poisson  equation.  The  particles  are 
"pushed"  in  the  calculated  electric  field  during  a  time  step  in  a  given 
cycle  and,  based  on  their  positions  at  the  end  of  the  time-step  new  charge 
densities  are  assigned  to  the  grid  points.  At  the  beginning  of  the  next 
cycle  a  new  potential  distribution  is  calculated  from  these  densities.  The 
calculation  thus  simulates  the  time-behavior  of  the  system  step  by  step. 

In  the  following  sections  we  describe  the  particle  injection  algorithms, 
trajectory  calculation  method,  and  technique  for  solving  the  Poisson  equa¬ 
tion.  Following  this  we  discuss  the  results  of  a  number  of  sample  runs. 

These  show  that  the  code  is  capable  of  handling  arbitrary  densities  and  tem¬ 
peratures  with  acceptable  fluctuations.  Current  collection  and  space  charge 
distributions  can  be  computed  as  functions  of  time.  Finally  the  program 
listing  of  PARKTDC  (Parker  Time  Dependent  Charging)  is  presented. 

The  present  program  is  an  original  development.  Although  there  is  a 
considerable  numerical  simulation  literature,  it  applies  to  other  geometries 
than  the  present  one.  Sources  available  to  the  author  give  little  or  no 
information  on  methods  of  simulating  axial ly-syn*netric  time-dependent  prob¬ 
lems  involving  finite  cylinders  in  space  plasmas.  Yet  the  axial ly-symmetric 
3-dimensional  geometry  requires  techniques  significantly  different  from 
those  of  say,  planar,  infinitely-long  cylinder,  or  cartesian  systems.  There¬ 
fore  no  references  are  cited. 
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PARTICLE  INJECTION  ALGORITHMS 

Particles  are  Injected  through  the  outer  grid  boundary  according 
the  following  velocity-distribution  options:  (a)  unidirectional  monoener- 
getic  beam,  (b)  isotropic  monoenergetic,  and  (c)  Isotropic  Maxwellian. 

Selection  from  isotropic  monoenergetic  velocity  distribution: 

The  outer  boundary  surface  of  the  grid  is  a  cylindrical  "box"  of 
finite  length  H  and  radius  RQ.  This  means  that  in  the  process  of  statis¬ 
tical  selection  of  particle  injection  sites,  one  of  the  three  areas  must 
be  selected,  namely  the  top  (A),  side  (B),  or  bottom  (C).  Consider  a 
distribution  of  particle  beams  distributed  in  polar  angle  0  (with  respect 
to  the  cylinder  axis),  at  a  fixed  energy.  For  an  isotropic  angle  distri¬ 
bution,  choose  e  from  cose  uniformly  distributed  in  the  range  (-1,  1). 
Define  the  projections  of  the  cylinder  areas  onto  a  plane  perpendicular 
to  the  beam.  There  are  two  cases: 

Case  e  >  tt/2 

Define  A  =  irRo^ | cose | and  B  =  2RQ  Hsine,  and  choose  random  number  to 
select  A  or  B. 

Case  e  <  tt/2 

2 

Define  C  =  irR„  cose  and  B  *  2R„  Hsine,  and  choose  random  number  to 
select  B  or  C. 

With  angle  e  and  areas  selected: 

If  the  selected  surface  is  A  or  C,  choose  azimuthal  angle  $  uniformly 
2  2 

from  (o,ir),  and  choose  r  uniformly  from  (0,RQ  ).  This  gives  the  position 
of  injection  on  the  surface.  If  the  surface  selected  is  B,  then  choose  z 
uniformly  from  (0,H),  and  choose  h  uniformly  from  (0 ,RQ) .  Then  obtain  ♦ 
from  arc  sin(h/RQ). 

It  should  be  noted  that  alternatively  one  could  have  first  selected 
the  areas  A  and  B  (or  A  and  C)  on  the  basis  of  a  uniform  distribution  of 
values  of  the  integral 

1  C  sine  1  cose |de 

tt  J  jcose)  +  (2h/irK0)sine 


(1) 
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which  represents  the  average  over  angles  e  of  the  ratio  of  A  to  A+B,  but 
which  does  not  seem  to  be  evaluable  in  analytical  terms.  Following  this 
selection  one  would  then  select  a  value  for  8.  This  method  would  be  equiva¬ 
lent  to  the  one  used. 

Components  of  velocity  at  Injection  point: 

For  a  given  speed  v,  and  for  selected  angles  e  and  <t>,  we  have  the 
components  of  particle  velocity  at  the  point  of  injection: 

x  =  v  sine  cos$ 

y  =  v  sine  si n<|>  (2) 

z  =  v  cose 


Selection  from  Maxwellian  velocity  distribution: 

Let  the  velocity  distribution  f(v)  be  described  by 

2  2 
3-c  -v2 

f (v)d  v-e  vxdvx  *  e  dvz 


(3) 


in  terms  of  dimensionless  velocity  components,  where  vz  is  the  axial  com- 

ponent  of  velocity,  and  v  is  the  perpendicular  component  /(i2  +  y2). 

"L'  p 

Then  choose  vA  from  RAN^  =  exp  (-Vj.  ),  and  vz  from  RAN2  =  erf(vz),  where 
RAN-|  denotes  a  random  number  uniformly  distributed  in  the  unit  interval, 
while  RAN?  denotes  a  random  number  uniformly  distributed  in  the  range 
(-1,1).  For  vz  we  need  the  inverse  (=erf  )  of  the  error  function  erf. 
(The  inverse  function  was  constructed  by  suitably  modifying  a  fitting 
formula  given  by  Abramowitz  and  Stegun.)  Then  the  angle  0  and  speed  v 
are  given  by  tana  =  vx/vz,  and  v  +  vz2.  Having  the  angle  e,  we 

then  select  the  injection  point  as  given  earlier,  and  following  this  the 
velocity  components  x,  y,  and  z.  If  there  is  also  a  drift  velocity 
M  (=  Mach  number)  in  the  axial  direction,  then  vz  can  be  chosen  from 
RAN 2  =  erf(vz-M).  This  method  of  selection  thus  leads  to  a  simplification 
when  dealing  with  a  drifting  Maxwellian  distribution. 


PARTICLES  INJECTED  IN  A  TIME  STEP 

The  actual  charge  of  particles  entering  a  surface  of  area  A  in  time 
At  is  ejQAAt,  where  j  is  the  random  thermal  flux  of  incident  particles 
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and  e  is  particle  charge.  Let  jQ  be  given  by  nQvA,  where  nQ  is  the  plasma 
density,  and  vx  is  the  mean  velocity  component  of  incident  real  particles 
normal  to  the  surface.  For  a  choice  of  time  interval  at =az/vL,  where  az 
is  a  chosen  thickness  of  transit  (say,  of  the  order  of  a  mesh  interval), 
we  may  determine  the  number  N  of  real  particles  injected  during  the  time 
step: 

N  =  jQAat  *  nQAaz  (4) 

If  n'  is  the  chosen  number  of  computer  particles  injected  in  a  time 
step,  say  10  to  1000,  the  charge  e'  per  computer  particle  is  given  by 


i 

e 


en.Aaz 
Ne  _  o 

W  N1 


(5) 


2 

In  the  problem  of  isotropic  Injection,  A  may  be  evaluated  as  2ttRq  (1+H/Rq). 
The  value  of  e'  is  used  to  calculate  the  charge  assigned  to  a  grid  point 
through  its  product  with  the  number  of  computer  particles  found  at  a  given 
instant  of  time  to  be  associated  with  that  grid  point.  It  is  similarly 
used  to  compute  the  charge  incident  on  and  absorbed  by  the  body. 


Typical  values  of  the  parameters  which  have  been  used  in  tests  are 
the  following: 

RQ  =  100  cm  N*  *  10  to  1000 

H  =  100  cm 


AZ  *  2  cm 

nQ  =  103/cm3  to  10^/cm3 
3  8 

For  nQ=10  ,  N=2. 51x10  is  the  number  of  real  particles  injected  per  time 
step.  Each  computer  particle  then  represents  N/N '  real  particles,  or 
2.51  x  105  to  2.51  x  107. 
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TRAJECTORY  CALCULATION 

The  particles  are  moved  in  accord  with  the  equations  of  motion,  which 
may  be  represented  as  follows. 

Let  X  denote  the  dimensional  vector  position  (X  stands  for  any  one  of 
the  3  cartesian  coordinates).  Let  T  denote  the  dimensional  time,  and  let 
V  denote  the  dimensional  potential  (in  volts)  at  the  position  X.  Then  over 
a  short  step  AT  with  constant  acceleration,  X  changes  in  accord  with 

x  =  Xo  +  xoAT  ‘  2m  (6) 


where  q  and  m  are  the  particle  charge  and  mass,  and  where  XQ  is  the  initial 
value  of  X  at  the  beginning  of  the  time  step.  Let  dV/dX  be  in  volts/cm,  and 
let  AT  be  expressed  as 


aT  = 


<aZ>min  at 
voi 


(7) 


where  (AZ)min  is  the  minimum  zone  thickness,  voi  is  the  scale  velocity  of 
the  ions,  and  At  is  a  dimensionless  time,  called  "DELTA"  in  the  program. 

Let  m  be  the  mass  of  the  particle  expressed  in  units  of  the  ion  mass.  Then 
we  may  write 


•  ^min  At 

X  =  X  +  X  - — - 

A  0  Ao  V. 


<AZU  <At>‘ 


01 


4  (m/m..) 


d 

dX 


where  E^ev)  is  the  scale  energy  of  the  ions  (=  m^v^- 
in  the  program,  and  qV  is  the  potential  energy  in  ev. 
the  velocity  components: 


gV' 

■V 


(8) 


72)  called  "TVIONS" 
Similarly  we  obtain 


X  = 


01 


^AZ^min  d  fqV' 

2 (m/m. )  31  |£7 


(9) 


Hence,  assuming  electrons  and  one  species  of  ion,  the  program 
(a)  reads  in: 

At  =  "DELTA" 

Ei(ev)  =  "TVIONS"  =  ion  energy 
Eg(ev)  =  "TVELEC"  =  electron  energy 
m^/me  =  "XMASS"  =  mass  ratio 
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and 

(b)  constructs  the  scale  velocity 

vQ  =  "SPEED"  *  5.93  x  107/TVI0NS/XMASS 


and 


=  ion  speed  in  cm/sec 


(c)  finds  the  minimum  zone  thickness  (aZ)m1n. 

The  particles  are  moved  in  subroutine  "TRACK"  (called  by  "DENSTY"),  with 
given 

X,  X/SPEED,  and  "DT"  »  (aZ)m1n  x  DELTA 

With  the  potential  grid  replaced  by  $a(mi/m)qV/Ei ,  the  new  values  of  X  and 
X  are  given  by: 


X  = 

X  * 


+  X 

-  v 


'  DT  (DT) 

o  v„  4 
0 

DT  <4 
o  T  dX 


2 


(10) 

(ID 


Interpolation 

The  required  components  of  d$/dX  in  the  latter  equations  are  obtained 
by  double  linear  interpolation  within  the  boxes  of  the  grid.  Let  r  and  z 
be  located  in  the  ranges  rjSr<rj+i  and  z^zsz^-j.  Then  the  interpolated 
values  of  a <j>/ar  and  a$/az  are  given  by: 

ff  3  [♦(U+l)  -  ♦(1J)  +  (z-z^)G/Dz]/Dr  (12) 

|f  3  frd-I.J)  -  4> ( "> » J )  +  (r_rj )G/Dr]Dz 

where 

G  -  ♦O-l.j+l)  +  ♦(l.j)  -  ♦(1-1J)  -  ♦(IJ+l) 
and 

Dr  =  rj+l  -  rj-  Dz  ■  *1-1  -  *1 

The  interpolated  potential  Itself  Is  given  by 

4>  3  *d»j)  +  (r-rj)[*(1,J+l)  -  ♦(1,j)]/Dp 
+  (z-z^OO-lJ)  -  ♦(1,J)3/DZ 
+  (r-rjHz^  )6/DrDz 


(13) 

(14) 

(15) 

(16) 


*V9' 
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THE  POISSON  PROBLEM:  POISSON  DIFFERENCE  EQUATIONS 


In  the  present  problem  the  electrostatic  field  is  axially  symmetric 
and  is  defined  on  a  mesh  of  spatial  grid  points,  such  that  at  any  point 
(including  grid  points)  the  potential  and  electric  field  can  be  obtained 
by  interpolation. 


Assume  that  the  space  charge  density  is  known  at  the  grid  points. 
Consider  a  group  of  interior  grid  points,  forming  a  portion  of  the  over¬ 
all  grid  as  shown  in  Fig.  3.  In  this  figure,  the  vertical  and  horizontal 
directions  are  the  z  and  r  directions,  respectively,  where  z  and  r  denote 
the  cylindrical  axial  and  cylindrical  radial  coordinates,  respectively. 


Three  horizontal  grid  lines,  of  constant  z-values  z._^,  z^,  and  z^,  and 
three  vertical  grid  lines,  of  constant  r  values  r^,  and  ,  are 
shown  in  the  figure.  (Note  that  the  index  (i)  of  z  increases  as  z  decreases.) 


The  set  of  grid  lines  intersect  at  9  grid  points,  or  nodes,  as  shown.  Each 
point  may  be  considered  to  be  associated  with  a  volume  of  space,  and  to  have 
a  group  of  four  neighboring  points  which  "interact"  with  it.  Thus,  consider 
the  central  point  of  the  group,  labeled  C  in  the  figure,  which  may  be  iden¬ 
tified  with  one  of  the  grid  points  in  Fig.  2.  Associated  with  this  point  is 
a  volume  of  revolution  (a  torus)  whose  cross-section  is  rectangular  and  is 
shown  by  the  rectangular  shaded  area  surrounding  Point  C.  The  shaded  area 
is  defined  by  connecting  the  mid-points  of  the  surrounding  mesh  rectangles. 
Let  t  denote  the  volume  of  the  torus,  and  let  the  neighboring  points  (above, 
below,  to  the  right  of,  and  to  the  left  of  C)  be  labeled  N,  S,  E  and  W 
(north,  south,  east  and  west,  respectively). 


Let  the  Poisson  equation  be  written 


(17) 


where  p  denotes  the  space  charge  density,  and  e  denotes  the  dielectric 
constant. 
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The  grid  lines  may  be  considered  to  be  arbitrarily  chosen  so  that  the 
mesh  intervals  are  nonuniform.  In  this  case  the  Poisson  difference  equa¬ 
tions  may  be  obtained  by  Integrating  Eq.  (17)  over  the  volume  r  of  the  torus 
associated  with  Point  C: 

J  |  dr  -  -  j  |  |  pdr  -  Qc/e  (18) 

T  T 

where  Qc  is  known  at  the  grid  point  C.  The  right-hand  side  has  been  approxi¬ 
mated  as  shown  since  x  is  small  In  principle,  and  Q  is  the  net  charge  in  the 
box  at  Point  C.  By  the  divergence  theorem,  the  left-hand  side  becomes 


where  z  denotes  the  surface  of  the  torus;  3<|>/an  is  the  component  of  v<j>  in 
the  outward  normal  direction  at  the  surface;  A^,  A$,  A^,  and  Ay  denote  the 
areas  of  the  north,  south,  east,  and  west  surfaces,  respectively;  and  the 
quantities  (a<j>/an)N  s  E  M  denote  values  of  a<|>/an  taken  to  be  constant  on 
the  corresponding  surfaces. 

(9<f>/an)N  s  E  w  may  be  approximated  by  difference  quotients,  namely. 


where  $  denotes  the  potential  at  Point  C  and  <1^,  denote  the 

neighboring  potentials.  If  Point  C  is  an  interior  grid  point,  the  areas 

AN*  AS’  AE*  and  AW  are  9iven  by 
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an  ■  f  C<r4+1  *  rj>2  •  *  rj-l>23 

AS  ’  *N 

AE  ’7<rJ*l  *  rj><z1-l  '  ZH1> 

Aw’f(rj*rM)(21-l-zHl) 

and  the  volume  t  Is  given  by 


T‘zi. 


'  zi*l> 


(21) 


(22) 


Thus  we  obtain  the  difference  equation  In  the  form 


CN  *N  +  CS^S  +  CE*E  +  *"  c4»  *  *”  Qc/6 


(23) 


where 


C  3  CN  +  CS  +  CE  +  CW 


(24) 


and 


_  MN 

r  -  __ 

N  U,_,  -  2,) 

cs  Tz. 

II 

m 

r  «s  - 

E  (rj+i  -  v 

CW  Jr 

1  ’  1+1 


(25) 


This  shows  how  to  form  the  difference  equations  used  for  the  Poisson  prob¬ 
lems  of  this  report.  Equation  (24)  holds  only  for  an  "interior"  point  of 
the  grid,  that  is,  a  point  surrounded  by  neighbors  on  all  four  sides.  If 
Point  C  has  a  known  neighboring  potential  (for  example,  if  Point  C  is  adja¬ 
cent  to  the  spacecraft  surface),  then  the  corresponding  term  on  the  left- 
hand  side  of  Eq.  (23)  is  transferred  to  the  right-hand  side  as  a  known  quan 
tity. 
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The  boundary  conditions  for  the  potentials  in  the  Poisson  problem  are 
as  follows.  At  points  representing  the  faody  surface,  the  normalized  poten¬ 
tials  are  fixed  at  the  chosen  values.  At  the  external  (boundary)  points  of 
the  grid,  where  "infinity"  is  represented  on  the  computer,  a  "floating"  con¬ 
dition  is  optionally  used,  namely,  a  linear  relation  between  $  and  3$/3n, 
the  normal  component  of  v$.  The  exact  relation  of  4  to  3<f>/an  is  not  impor¬ 
tant  when  the  external  boundary  of  the  grid  is  sufficiently  far  away.  (For 
the  calculations  to  be  reported,  the  assumed  relation  was  the  same  as  for  a 
Coulomb  potential.)  In  any  case,  either  the  fixed  condition  $  =  0  or  the 
floating  condition  will  give  the  same  results,  provided  the  grid  boundary 
is  moved  sufficiently  far  out.  The  effects  of  various  types  of  boundary 
conditions  representing  "infinity"  have  been  studied  by  the  author. 

In  general,  the  floating  condition  appears  to  be 
computationally  more  efficient  than  the  fixed  one.  Of  course,  the  floating 
condition  becomes  ideal  when  the  true  relation  between  $  and  3<|>/3n  is  used, 
but  this  requires  that  the  asymptotic  form  of  the  solution  be  known  in 
advance.  The  boundary  conditions  at  the  outer 

grid  surfaces  can  be  combinations  of  fixed  and  floating  conditions. 

Consider  a  Point  C  on  the  outer  boundary  of  the  grid  where  a  floating 
boundary  condition  is  chosen.  If  the  potential  is  assumed  to  satisfy  the 
linear  law 

3<t>  _  d#  _  /oc\ 

-  -§z  “  “  “♦  (26) 

on  the  z-boundary  (North  or  South),  and 

|i  =  |i=-W  (27) 

on  the  r-boundary  (East  only;  6*0  on  the  West),  then  the  correspondi ng 
"neighbor  term"  on  the  left-hand  side  of  Eq.  (23)  vanishes,  and  the  cor¬ 
responding  "neighbor  coefficient"  on  the  right-hand  side  of  Eq.  (24)  is 
replaced  by  aA  or  SA,  where  A  Is  the  appropriate  area.  The  quantities  <* 
and  0  depend  on  the  position  and  on  the  assumed  model  for  the  variation 
of  the  potential  at  large  distances. 


Once  the  coefficients  of  all  of  the  equations  (corresponding  to  the 
grid  points  where  the  potentials  are  unknown)  are  computed,  the  system  of 
linear  equations  of  the  form  of  Eq.  (23)  may  be  solved  by  iteration.  Point- 
successive  over-relaxation  is  a  well-known  process  and  has  been  found  to  be 
effective  in  the  present  problem.  For  the  relaxation  process,  one  rearranges 
the  equations,  so  that  the  "diagonal'1  term  is  alone  on  the  left-hand  side, 
while  all  the  other  terms  are  on  the  right-hand  side  with  the  known  charge- 
density  term.  Thus,  Eq.  (23)  becomes 

C*  ■  ^  *  Qc/e  (28) 

First,  an  initial  guess  is  made  for  the  values  of  all  the  potentials.  Then 
new  values  are  obtained  from  the  left-hand  sides  of  all  of  the  equations 
(28),  using  previous  values  on  the  right-hand  sides.  One  "sweeps"  through 
the  equations  successively,  replacing  the  potentials  on  the  right-hand  sides 
with  updated  values  as  they  become  available  from  preceding  equations.  This 
procedure  is  usually  stable  and  leads  to  convergence.  "Over-relaxation"  is 
the  process  of  mixing  successive  potential  iterates  in  such  a  way  as  to 
enhance  the  rate  of  convergence. 

It  is  convenient  to  express  all  potentials  in  volts  and  all  lengths 
in  centimeters.  Then  if  the  charges  are  all  expressed  in  picocoulombs,  we 
need  simply  to  replace  1/e  by  0.9x4*  or  3.6*,  in  vacuum. 
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SAMPLE  RESULTS 

One  of  the  most  vexing  problems  of  numerical  time  simulation  is  that 

of  noise,  wherein  fluctuations  not  representative  of  the  actual  plasma  are 

generated  by  the  relatively  small  numbers  of  particles  and  grid  points  used. 

This  problem  becomes  more  severe  as  the  plasma  density  increases,  that  is, 

as  the  Debye  length  becomes  small  compared  with  the  spacecraft  size.  Since 

it  was  deemed  essential  that  the  code  be  able  to  cope  with  dense  as  well  as 

tenuous  plasmas,  as  part  of  the  code  checkout  sample  runs  were  made  with 

3  3 

the  parameters,  temperature  T*0.1  ev,  and  density  nQ  varying  from  10  /cm 
to  as  high  a  limit  as  practicable. 

Among  the  quantities  of  interest  are  the  collected  fluxes  of  ions  and 
electrons.  Fluctuations  in  these  are  a  measure  of  those  in  the  plasma.  In 
the  first  problem  to  be  discussed  the  spacecraft  is  modeled  by  a  disk,  of 
diameter  one  meter.  The  plasma  has  density  1000/cm  ,  and  is  Maxwellian  with 
T=0.1  ev.  In  this  case  the  Debye  length  is  7.43  cm,  so  that  the  Debye  num¬ 
ber  is  0.149.  The  electron  and  ion  plasma  periods  are  3.5xl0"6  sec  and 
1.8xl0"5  sec,  respectively,  for  an  assumed  ion  mass  of  25  electron  masses. 

(The  use  of  an  unphysically  small  ion  mass  is  common  in  simulations,  espe¬ 
cially  if  the  steady  state  is  wanted.  In  the  steady  state  the  solution 
normally  does  not  depend  on  the  ion  mass.) 

The  disk  is  assumed  to  "turn  on"  at  t=0  with  a  fixed  potential  of  -10  v. 
Time  steps  of  length  2. 36x1 0"6  sec  are  chosen,  with  100  electrons  and  20  ions 
injected  per  step.  (The  ratio  100  to  20  is  the  same  as  the  ratio  of  electron- 
to-ion  random  thermal  fluxes.)  A  coarse  grid  of  25  points  is  used. 

The  number  of  ions  absorbed  per  step  represents  the  collected  ion  cur¬ 
rent.  This  quantity  is  shown  in  Fig.  4  as  a  function  of  cycle  (step)  number. 
It  is  strongly  fluctuating,  but  quite  clearly  centered  about  12/step. 
Increasing  the  number  injected  per  step  to  1000  electrons  and  200  ions 
results  in  a  smaller  fluctuation  as  illustrated  in  Fig.  5.  When  10-cycle 
averages  are  plotted  instead  of  the  raw  data,  one  obtains  a  "quieter"  current 
as  shown  in  Fig.  6. 

4  3 

In  further  runs  to  be  discussed,  the  plasma  density  varies  from  10  /cm 
to  105/cm3.  The  corresponding  values  of  Debye  length,  Debye  number,  electron 
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plasma  period  xe,  and  the  ion  plasma  period  ,  are  given  in  the  following 
table. 


Run 

(cm“3) 

,XD 

(an) 

X0/50 

Te 

ik 

ill 

A 

104 

2.35 

.047 

1.11  X 

10’6 

5.55  x  10‘6 

B 

2.5  x  104 

1.49 

.030 

7.02  x 

io-7 

3.51  x  10‘6 

C 

5.0  x  104 

1.05 

.021 

4.96  x 

10"7 

2.48  x  10-6 

D 

7.5  x  104 

0.86 

.017 

4.05  x 

10'7 

2.03  x  10-6 

E 

105 

0.74 

.015 

3.51  x 

10-7 

1.76  x  10"6 

In  the  following  we  will  discuss  Runs  A  and  D  in  detail.  Runs  B  and  C 

gave  results  intermediate.  Run  E  was  so  strongly  fluctuating  that  it  was 

4 

not  pursued  further.  The  density  of  7.5  x  10  seemed  to  be  the  largest 
(at  T=0.1  ev)  that  could  be  treated  at  the  present  boundary  condition  (20  cm 
from  the  spacecraft  surface).  Higher  densities  will  require  that  the  boun¬ 
dary  be  moved  inward. 

Figure  7  (Figs.  7a-7g)  shows  results  of  Run  A  (n  =104).  This  is  a 

0  -8 

printer  plot  of  the  variations  with  cycle  number  (one  cycle  =  9.5  x  10  sec) 
of 

A  =  ion  absorptions  per  step  (scaled  from  0  to  30) 

B  =  ion  population  (scaled  from  2500  to  5000) 

C  =  electron  population  (scaled  from  2500  to  5000) 

D  =  maximum  potential  (scaled  from  0  volts  to  25  volts) 

The  significances  of  A,  B,  C,  D,  are  as  follows.  We  inject  100  electrons 
and  20  ions  per  step. 

A  represents,  through  the  ratio  of  the  number  absorbed  to  the  number 
injected,  multiplied  by  the  ratio  of  the  area  of  collection  to  the  area  of 
injection,  the  collected  current  relative  to  the  ambient  thermal  current 
available. 

B  and  C  represent  the  ion  and  electron  populations  (numbers  of  parti¬ 
cles  active  throughout  the  grid).  These  should  become  constant  as  the  steady 
state  is  approached.  In  order  to  arrive  at  the  steady  state  as  quickly  as 
possible,  a  "quiet  start"  approach  is  used,  wherein  the  space  is  previously 


FIG.  7a. 


FIG.  7d 


filled  with  particles  so  as  to  represent  normal  density.  In  the  present 
case  5000  ions  and  5000  electrons  are  used,  to  represent  the  unperturbed 
populations. 

0  is  a  measure  of  the  fluctuation  intensity.  Since  the  boundary  poten¬ 
tials  are  allowed  to  float,  they  are  susceptible  to  oscillations.  In  par¬ 
ticular,  when  the  instability  is  severe  the  boundary  potentials  tend  to 
become  positive,  with  the  maximum  Increasing  as  the  severity  increases. 

Hence  the  maximum  is  monitored. 

The  behavior  of  A,  B,  C,  and  D  as  functions  of  time  are  shown  by  Fig.  7 
as  follows.  Large  transients  occur  in  all  4  quantities  during  the  first 
150  cycles  (t=1.4xl0‘5  sec).  Both  the  electron  (C)  and  ion  (B)  populations 
decrease  from  their  unperturbed  values,  the  electrons  because  they  are 
repelled  away,  and  the  ions  because  they  are  accelerated  and  spend  less  time 
in  any  given  volume  (a  consequence  of  Langmuir  probe  theory).  The  electrons 
(C)  oscillate  while  approaching  their  final  value,  the  oscillations  having 
diminishing  amplitude.  The  period  of  the  oscillations  is  roughly  50  cycles, 
or  4.7x10"®  sec,  about  4  times  the  electron  plasma  period  (see  table), 
approximately  equal  to  the  ion  plasma  period.  It  is  not  clear  whether  the 
oscillations  represent  physical  behavior;  their  smoothness  suggests  that 
they  do.  However,  the  effects  of  numerical  "aliasing",  wherein  high- 
frequency  sources  can  drive  low-frequency  oscillations  through  the  finite 
grid  spacing,  must  certainly  be  present  because  of  the  relatively  large  grid 
intervals  (of  order  10  cm). 

After  about  220  cycles  (2.1xl0“5  sec)  the  electron  population  has  set¬ 
tled  down  and  subsequently  rises  slowly  toward  its  asymptotic  value,  about 
3500  out  of  the  original  5000.  The  ions  (B),  on  the  other  hand,  do  not  oscil¬ 
late,  but  drop  monotonically  toward  a  shallow  minimum  at  about  180  cycles 
(1.7xl0-5  sec),  and  then  rise  slowly  toward  their  asymptotic  population, 
about  4400  out  of  the  original  5000. 

The  ion  absorption  rate  (A)  rises  Initially,  overshooting  and  subsequently 
decaying  to  an  asymptotic  value  of  about  11  absorbed  per  step.  A  higher- 
frequency  oscillation  of  relatively  small  amplitude  and  a  period  of  about  20 
cycles  (1.9x10"®  sec)  is  superimposed.  The  ratio  of  the  current  collected  to 


29 


the  unperturbed  current  that  would  be  collected  In  the  absence  of  electric 
fields  is  given  by  dividing  the  number  absorbed  per  step  by  the  number 
injected  per  step,  multiplied  by  the  ratio  of  injection  area  to  collection 
area,  and  divided  by  DELTA,  a  time-step  input  defined  earlier. 

Figure  8  (Figs.  8a-8g)  shows  results  of  Run  D  (nQ*7.5xl04) .  All  inputs 
are  the  same  as  in  Fig.  7.  However,  now  we  see  that  there  are  severe  fluc¬ 
tuations  in  curve  D,  the  boundary  potential.  The  ion  and  electron  popula¬ 
tions  (B  and  C)  both  descend  almost  monotonlcally  to  lower  values  than  they 
had  in  Fig.  7.  They  are  closer  to  each  other  (about  3300  and  3000)  than 
they  were.  The  ion  absorption  rate  (A)  initially  overshoots  to  a  large 
value,  but  more  quickly  settles  to  a  ''steady1'  state,  with  some  fluctuations 
though  not  severe,  about  6  absorbed  per  step,  roughly  half  of  that  in  Fig.  7. 

During  the  transient  period  in  absorptions  (first  70  cycles  or  about 
6.6x10"®  sec),  the  electron  population  (C)  oscillates  with  relatively  small 
amplitude  and  smaller  period  (about  25  cycles  or  2.4x10"®  sec)  than  in  Fig.  7. 
This  period  is  again  about  the  same  as  the  ion  plasma  period  (table).  Hence 
this  oscillation,  which  appears  regular,  may  be  physical.  It  is  interesting 
that  the  electron  population  fluctuations  (the  ions  have  little  or  none)  are 
much  less  severe  than  and  apparently  uncorrelated  with  the  boundary  potential 
fluctuations  (D). 

It  should  be  stated  that  a  "cloud-in-cell"  model  was  also  developed  and 
used,  without  significantly  reducing  the  noise  as  had  been  hoped.  It  also 
gave  slightly  different  (higher)  values  for  the  steady  state  absorption  rate. 
The  runs  made  with  the  cloud  model  will  not  be  discussed. 

Figure  9  shows  a  potential  contour  plot  for  Run  D,  generated  by  the  com¬ 
puter  printer,  at  the  "steady"  state. 
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PARKTDC  PROGRAM 


The  listing 
pages.  (PARKTDC 


for  the  PARKTDC  computer  program  is  given  in  the  following 
*  Parker  Time  Dependent  Charging.) 
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— ARE  wtif  J)  »r  I  *t»c  2**c'  " — f - - - 

NCH6€tJ)=0 

- - 

171 

IF(IT.EQ.O.AMO.IRELAX.6T.O)  PHICJ>=0. 

- t/g - 

173 

c 

175 

Ufnimil  •  ’  *  3  '  ■  . . 

00  525  1=1, N2A 

I  r  C 

177 

ZZ— Zfllll)  - 

IF(I*LT#NZA)  Z1-7AHI+1) 

tr  e 

179 

0EL2CI)=  22  -  Z1 

100 

161 

525 

CONTINUE 

. .  iO  c 

163 

C 

XHASS=HASS 

to* 

16  5 

Ii  . 

00  550  I«1,NZ1 

187 

Icc*i££^l  . T  1  "  -  ..  —  ■ 

Z2=Z0I€I1 

XO  O 

189 

------  4  on 

IrTXtt **N£OI 

IFd.EC.NZ3)  Z1=ZBICIZBJ 

l  yO 

191 

ecu  _ 

'OLtCiiCCF*  is.  cl 

0ZMIN=AMIN1 ttZMIN, OELZ  (IZZ)  ) 

193 

1  44 

c 

TPTDfTafl  - 

195 

4  Qf 

_  ..  i* 

lr  tH 3f*f  .  .  - 

JFIRST=0 

197 

_ ±J1A _ 

c 

rtONHAX=0  F Ofi  MONOENERGETIC  BEAM 

iiAuMiv.t  ene  reflTBABTP  urturtCkiforcTTA 

199 

_ 243 _ 

t 

MGNH AX=2  FOR  ISOTROPIC  HAXMELL  IAN 

201 

_  4(1  p 

y 

H0NMAX=2 

fcllfc 

203 

IF  t  H0NHAX.EQ.1)  0EN0M=4. 

20  5 

c 

20  T 

_ 3ft  A  _ 

IFtHONHAX.GT.Or  PCON=PCONM1.*  ill  Cl) -ZZ INZI)  /RR  CNP  )  )*2.  209 

_ ■* _ gf-TJftAAkl _ HJl _ 

. 

XIINJ»NIINJ 

KnuBiB.apnuivTTu  r 

211 

_ 2-12 _ 

SPEE0E=5 • 93£7*SQRT (TV ELEC) 

213 

°ROGRAM  PARKTDC 


12/11/00  PAGE  5 


SPE£DI=5 .93E7*SQ  FT (TVICNS/XMASS) 


iniili:,  t*  ;  i  >1 -a  -»  j  i  ■  r*ki  a;  .  i 


CURRE=1.  6E-7*DENCC*SPEE0I/0EN0M 


PL ASP£=l«llE-4/SQRT  CDENCC) 


OTSEC*OELTA*OZMIN/SPEEO*OENOM  224 


HRITEfM,  9920)  DELTA, OTSEC, TIME,  TVIGNS,TVELEa,OENCC , PLASPI ,PLA SPE, 


uyv  i  fcf r 

IF<MOOE.EQ.ll  MRITE(M, 9921)  WCLOUO  228 


NFPP= (NT0T/3G0) ♦  1  230 


WRITEtM, 90301  232 


60  0  CONTINUE 


C  INITIALIZE  MATRICES  FOR  ANALYSIS  OF  VARIANCES 


00  610  1=1, NTOT 


AVGE<I)»0# 


RMSOECI) =0* 


620  CONTINUE 


IFCIT.GT.01  GO  TO  650 


INITIALIZE  SPACE  CHARGE  MATRIX 


V0L=PI*RR«NR)**2*<ZZm-ZZ(NZn  2( 


MNJ=RPART/CCMPAR  2J 


w  wviinw  w  *  *  r 

WRITE (N, * )  "  INITIAL  Nl NJ,RPART, COHPAR  =  ", NINJ,RPART,C0MPAR 


5 


45 


program  parktoc 

12/11/8  0 

PAGE  6 

CALL  DISCON(M) 

265 

00  uJu  c 

IPART=IOUM 

- 2-6r€ - 

267 

"  “CuO*3PtcOI  .  ” 

IP C IPART • £0*2)  SPEEOO*SPEEO£ 

- 2t« 

269 

• 

IFI*ST=IFIRST*1 

- 2  79- . 

271 

OJv 

INIT=0 

272 

273 

4 

c 

65  0 

CONTINUE 

275 

c 

- B - 

COMPOTE  IPRINT 

- ere - 

277 

- - 

IPRINT=0 

279 

IF(NFRINT  .EQ*1)  IPkINT**! 

IF CHOOfIT $ JSKIP1 #LE«I SKI P-1)  IPRINT*! 

- 280 - 

281 

- fi - 

IFtIT.EQ.O)  IPRINT  =  3 

'  Z'UC 

283 

- — _ 

C 

* 

C 


IF(JFIRST.GT.OI  TIME* TINE+OTSEC  285 

CALL- ■■F-gfcfl - 29E 

IF < IPRINT aGT«l)  WRITE  CM* 9790)  IT  287 

EO> *)•  C*t4  CONNE04M) - 289 

CALL  FPL  OTCXfRRjZZ»NR»  NZ»  NRC*NZO»  RADIUS,  WAX*  IT)  289 

- 2-99 

IFUT.GE.ITHAX)  GO  TO  1000  291 

- : - - . - - - rr-r— • „ - 299 

293 


IPART  fcOOP - IPA6T»*,g,  A NO  8  FOR  IONS,  ELECTRONS?  AND  PHOTOELEGTR - 


295 

00  900  IPWlyfr---" - - « - 296- 

IPART=IOUM  297 

IFC  IPART,  Sq»i»T«OtTS»  THIONS — - - - - 299 

IF4IPART.EQ.2)  TVOLTS=-TVI0NS/XHASS  299 

IF41PABT w  EQ,1)  SP€€g0«a»€E01 - 399 

IFIIPART.EQ.2)  SPEEOO=SPEEOE  301 

IF  4  IPART.  E&w&r  -TWCHG»TWICNS - - - - - 392- 

IF  < IPART • £Q*  2)  TYCHG=-TVEL£C  30  3 

IFtIPART.EQ.l)  CURRiCURRI - 300 

IFf IPART •  EQ« 2)  CURR=CUfiRE  30! 

IF4  IPART.  EQ.l)  P ARTCL <  1  )~»P  ART1 449 - 398 

IF(IPART-EQ.l)  PARTCl  «2>=PART1<2>  307 

IF«IPABTttGT»1)  PARTOL  4-i>»PART2f*) - 309 

IFIIPART.GT.1)  P ARTCL  C 2)=PART2 12)  309 


• 

- 6 - - — • — - — — — - 

r,***+  + 

NEXT  TIME  *T£E_l 

-3*0 - 

311 

C  -  -------- 

313 

IF 4 IPRINT •GT»1)  WRITE 4M ,9750)  PARTCL, SPEEDO 

yCUKR 

315 

6 


PROGRAM  PARKTOC 
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C 

- WIMJaWTINJ - ! - - - - 

IF ( IF ART •  EQ*  2)  NINJ«XIIN4*SPEEOE/SPE£OH>.1 

tj - ( - - - : - - 

CALL  OENSTY(INIT) 

t - ; - 

C  COAHB=2«  677*1TENCC*SQRT <T  VELEC1 

■c - - - - - 

J¥AR=MOO t JFIRST /2» 10) *1 

- PO  TOO  N-l, WTOT - — 

IFdPART  .EQ.2)  GO  TO  690 

- 0  ( N  )-PKOUNt  N)  »CHeO N/XITN  J - : - 

NETCH(N) =KOUNT IN ) 

- KHISTTtttf  tl»-KHI3TI(N ,  11) ♦KOUNTtN) KHttttm 39**1 - 

KHISTI  (N » 12 )=KHISTI (N » 12) +KOUNT( N) *KOUNT (N)  - 

- 1  KHlSTItm-JVARt  »KHI  ST  I  ( W  ,  JVAR) - ■ - 

KHISTI (N , JVAR)=KOUNT(N) 

- 00  TO  TOO - - - — - 

690  CONTINUE 

- o  tn )  -o  cm  ■rKcuntfo  *CHCtsnttftm - 

NETCH  (N) =NETCH( N) -KOUNT(N) 

- ICHI?TetNill)-KHISTgm»ltr’*'KOUNTtHr’»imtSTgtNi  JtTAItl - 

KHISTE (N»12)=KHISTECN»12)*K0UNT( Nl *KOUNT (N> -KHISTE ( N, JVAR) * 

- 1  KHISTgCWiPPAgl - - — - - - - - = - - - - - 

KHISTEIN ,JVAR) *KCUNT(N) 

— Ypfl - CONTItrtJg . . = - - s - — - . . . 

C 

P - APP  TO  REMAINING  VARIANCE  SUMP - — - - - — - 

C 

- If (I»ART.C&.g)  OP  TO-  »gO — - - - - - 

NIA8S=NCA0S 

- NIACT^NPAC-T . .  . . - — — - 

NIESC=NCESC 

- KHA03miT=KHA9SI(ll>»NCAe3’-XHA93I(  JVAR) - 

KHACTI (lllsKFACTI (11) +NCACT-KHACTI ( J¥ AR) 

- KHAtH»l(tg»^«mAaPI-<ie»*NeAP9«HeAa9-KHA93I(JWAR)>KHAPSI<JVAR) 

KHACTI(1Z)=KHACTI  (12)  ♦NCACT*NCACT-KHACTI(JV AR)  *KHACTI  ( JVAR) 

- KHAg3I(JVAR)*NCA93 - - - — — - - - 

KHACTI ( JVAR)=NCACT 

- GO  TO  798 - — — — — - - — - 

720  CONTINUE 

- WgAP3-«NPA93 - 

N€ACT=NC ACT 

- NggSC^NCgSC - - - - - 

KHA6SEI11)=KHABS£( 11) ♦NCAgS-KHABSE ( JV AR) 

- KH A CTg(t~lP*  RH ACT g <  1 1)  ♦NCAeT-RHACTgCJPR*) - 

KHA8SE(12) =KFABSE( 12) ♦NCA9S*NCA8S-KHABS£< JVAR) *KHA9SE (JVAR) 

- KHAeT€(lg7^'KHAeTetlg)»NCACT»NCAeT-KHAeTgtJ«AR)  * KHACTC f"JPAR)~ 

KHAeS£(JVAR)sNCABS 

- KHACTgtJVA»)=NCACT - - 

72 0  CONTINUE 


3ie 

srr 

318 

STT 

320 

J9tf 

322 

P2P 

324 

1?t 

326 

P2P 

328 

*2* 

330 

psr 

332 

PSP 

334 

PSP 

336 

PST 

PSP 

340 

P4T 

342 

P4P 

344 

P4P 

346 

PAP 

348 

PAP 

350 


PSP 

354 


356 

PSP 

358 

PSP 

360 

per 


364 

PP5 

366 


JFIRST=JFIRSm 


NCHGSE=0 


00  71Q  J=l,  JEOGE 


NCHGSI=NCHGSI+NCHGI  <J) 


3 

— 3 
375 
--J76- 
377 


00  910  J=i, JEOGE 


PM I  RELATION  HILL  9E  ESTABLISHED  LATER*  THIS  IS  INCOSPECT. 


PHI (J)=PI/2.*.9*Q0ISK/RA0IUS 


920  CONTINUE 


00  930  1=1, NZ 


N2=N2*NR 


CONTINUE 


C  BEGIN  CALCULATION  OF  VARIANCES 

_e - . — - — - . - — - 

IF ( ISH6»  EQ* 1«  ANO* JFIRST*  EG*  2)  CALL  CONNEC(M) 


AVI  A9S=FL0  AT  (KHABSHil)  l/CELT 


RMIA3S=SQRT<FL0AT<KHABSI 412) l/OELT- AVIA3S*AV 3A9S) 


AVEABS=FL OAT (KHABSEf 11) l/OELT 


RMEA8S=SQRT  <FL0AT<KHABSE<121 1/ CELT- A¥EA3S*AV£ABS) 


IFCIPRINT *GT*0)  WRITE  CH,  3002) 


3002  FORMAT (///25X,21HSUHM ARY  FOR  THIS  STEP, 16X, 19HTEN-STEP  RMS  ERRORS, 

- 1  17V,  1 7HT E N—STEP  AVER  AGES/  22X ,  »HK)NSr-^X  ,BHEL£GTRONS ,  IX  ,  2  <  3X, 

2  1HI, 7X, 1HE) )/lX, 5HCYCLE , 3X,4HTIME*2f4X» 14H ASS  ESC  ACT), 


IPLAST=IPRINT 


1  RNIA3S,RMEA9S ,RMIACT , RNEACT, A VIA9S, A VEA3S, A VI ACT, AVEACT , XMAX 


PROGRAM  PARKTOC 
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WRITE(KV,3001>  IT,TIME,NIABS,NIE SC*NIACT,NEA8S*NE£SC,NEACTt 


30  01  FORMATtrX,I4,iP£18.3,iX,I3f'2I6,3X,I3,2I6,3X,OP4F8.1,3X,4Fa.i»F8.3> 


CO 


GO  TO  650 


99  CONTINUE 


*  f  ■  n  m  *  ■  K3  •  < 


REAOI  JV)  INJA, (tVSAVE <I»  J) »I=1, 10 > , J=l, I NJA) 


1030  CONTINUE 


GO  TO  10  20 


1010  CONTINUE 


CREATE  RECORC  TO  BE  COPIEO  TO  INPUT  FILE  FOR  RESTART 


REMINO  JV 


MRITEf JV,1111)  JEOGEtNR* NZA,NZ8, IRELA X,MOOE 


WRITE <JV*  2222)  <ZA  (I)  ,I=1,NZAI 


PROGRAM  oahxTOC 
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I 


WRITECJV, 2222)  (Gtl) , I=l,NTOT) 


STOFll 


FORMAT (1 615 ) 


3333  FORMATI1X,020) 


5555  FORMAT tiHi ) 


1X,I3,25H  COLUMNS  IR-VALUES)  TOTAL  /  487 


3  1X,I3,46H  ROWS  (Z-VALUES)  BELOW  ANO  INCLUDING  Z  =  0  PLANE  / 


6  10H  IS  FlXEC,,48H  *1  IF  800Y  POTENTIAL  IS  FLOATING,  ANC  EQUILI3R 


99 SI  FORMAT C/ IX, 13, 41 H  *  MOOE*  SELECTS  NORMAL  CALCULATION  MODE) 


9993  F CRM AT (/IX, 13, 41 ►  =  MOOE,  SELECTS  WEIGHTED  PARTICAL  HOOEL) 


1  I3,1PE15.4,3H  CM))  497 


FORMAT (// 34H  Z-VALUES  POSITIVE  A30VE  Z=0  PL ANE/ { IX , 13, 1PE15 . 4  , 

1  3H-  GH->) - — -  - - - 5*0- - 

9950  FORMAT (/✓ 34H  Z-VALUES  NEGATIVE  BELOW  Z=0  PL ANE/ ( IX , 13 , 1PE15 . 4 , 

- 1  3H  CM)) - - - - - = - 5-02 - 

9940  FCRMATI//2TH  SURFACE  POTENTIALS  ON  BODY/ CiX  ,  13,  lt>E15.  4  ,  503 

a  e  u  iimi  Te  %  a  _  _____  ... _ _ _ _ _ _ _ _ CiLL.  — 


9930  FORMAT  C//38H  INTEGER  INPUTS  FOR  SHEATH  CALCULATION  /IX,  505 

- -a— l3y37H-—  N PRINT  (TO  PRINT-  TRAJECTORY  STEPS Wi*, - 5<H - 

2  I3,39H  =  NPTS  CTO  COMPUTE  ONE  OR  MORE  POINTS)/lX,  50? 

- 4  2I5,44H  ■  IT  ANO  ITS  -< ITERATiQH--IH£>€X  -ANO-NUM^ER-OF, - 5*2 - 

5  12H  ITERATIONS) /IX,  509 


7  I5,35H  PAGES  OF  SUMMARY  PRINTED  FOR.  EACH, 15, 6H  STE^S,  511 

- a  13 , 3BH  »  MASS-4MA5S  OF  COMPUTER  -I 3W  ) - 51-2—— 

9920  FORMAT <// 27H  EMISSION  ANC  PLASMA  INPUTS  /  513 

- 1  lX,P10.3,3ftM»-OEbTA - -  INITIAL  TIME  STEP - CME^fS, - 514 - 

2  iPE15.4,3X,«HSEC0N0S> /  515 

^  x  •  .  «  «  j  ■  ■  «  «>.«•  _  a**  i  i  w  «  I 


4  1X,F10. 3,49H  =  TVIONS  *  AMBIENT-ION  TEMPERATURE  IN  VOLTS/ 

5  lX-,f-l0«  3, 4  OH  ~  TVELEC  t -AM9I£HT^6LgS-TRON  TEMPER  AT  URE— IN-  VGt-T-S/ - 

6  1X,F1Q.3,5C-  =  DENCC  *  AMBIENT-ELECTRON  DENSITY  IN  NO,  PER  CC  / 


10 


7  21)X«  25H  ( ION  PLASMA  PERICO  =  ,  1PE 15 • 4, 3X , 8F SECONDS) /  520 

- 8-  20XV25H  (ELECTRON  PLASMA  "PERIOD  a,E15,4,3X,  8HSSCONCS1/ - 5« - 

9  20 x, 2 5hc electron  oebve  length  =,eis.4, ex, 3hcmi/  522 

- a-lTttP-ErSTAiTOW  -v  PCON - »"NUHOER  ~0'F"RCA L  IONS  INJECTEO  PER  CYCLE/ 

a  IX, £15.4,53*  =  COMPAR  *  NUMBER  OF  REAL  IONS/ELECTRONS  PER  COMPUTE 

- C— iftfflT- PARTICLE) - m — 

9921  FCRNAT(1X,1P£15.4,29H  *  4CL0UD  =  CLOUO  HALF  WIDTH)  526 

"988D FORMAT  C/7~1'X ; 29H  INPUTS  FOR  SINGLE  TRAJECTORY/ - 5Tf - 

1  1X,F10. 3,19*  *  RPT  *  R-POSITION/iX,F10.3f19H  s  2PT  =  Z-POSITICN/ 

- 1  1X»F10»  Tj  3 OK  =  AL1  =  POL~AR"  ANGLE  (OgGREES)/ - 5T9 - 

3  1X,F10.3,34H  =  3E1  *  AZIMUTHAL  ANGLE  (DEGREES)/  530 

- *r  -  ~ex>  a  energy  in  volts) - «t — 

9870  FORMAT (//13H  ALL  Z-VALUES/ (IX, 13, 1PE15.4))  532 

^08U - FORMAT f/'/ IX  jtg'l HINT ERST IT IA~L  REVALUES/  (IX, 13, lPEtS.  4D - SS3 - 

9850  FORMAT  (//IX, 30HINTERSTITIAL  Z-VALUES  POSITIVE/  (IX,  13, lt>E  15.4)  ) 

■984TO - FORMAT  (Tiny  fTOHTNTERSTTTIAL'  Z-VAL  tlES  NEGATIVE/  tlX,T  SytPgtr.'  4  » - 

98 35  FCRMAT(//1X,43HINTERSTITIAL  Z-VALUES  NEGATIVE  ANC  POSITIVE/ 


9830 

1  (lX,I3^1af.l5.4>) - 

FORMAT (1H1 76X, 1HN, 5X, 4HR (N), 4X, 4HZ (N) //) 

- 9Tf - 

53  fl 

9T90 
978  0 

F  CRM  AT  (AIK  ?  X  3  $  16  M  CROCK  POTENTI A-L# 

F0RMAT(/1X,  13, 25H  OROER  FLUXES  ANO  CHARGES) 

539 

540 

975  0 

FGRRAT</i1tf  23HGOING  TMTO  uCMSTY  WITH  f  2A5f  5Xf 

1  13H  WITH  SPEEO  =,1PE15.4,7H  CM/SEC, 5X,22H  ANO  CURRENT 

541 

OENSITY  =, 

9740 

2  E15.4,14H  PICOAHP/CN* *2) 

FORMAT (//5X27HSUMMART  OF  NET  SPACE  CHARGE// 12X ,9 (A2 , 

543 

544 

97  38 

1 , 2HR (, 11 , 1H ) ) , 1 1 ( A 1,2 HR { , I Z »1H) ) )  - 

FORMAT (/5X,2*Z(,I2,1H) ,2016) 

. . .  545 

546 

51 


PROGRAM  PARKTOC  /LIST 


SUBROUTINE  LIST(LST,IP) 


PRINT  ARRAYS 


0 2/11/6  0  PAGE  12 


CCMNON  H,IV,IFIRST, JFIRST, JEOGE,NR,NZ A, NZ, NTOT.PI , IT , IR, 12, 


2  AREAHC1 1 ) , OEL Z(11),RI(11) , ZCI (11) 


00  200  NP=1, 5 


IFCKP.GT.NTOT.ANO.NP.EG.  1)  RETURN 


NMAX=NP 


IKP=(KP**1 )  /NR+1 


IFCLST.EQ.l)  ZOUTINP) =  X  t  KP) 


IF (LST«£Q«  2 J  20UT (NP) -22 IJKP) 


200 

_ _ rL_ 

CONTINUE 

572 

"300 

CONTINUE 

7  /  J 

574 

IT  1LJMI fCOll  - IV  WZT9 = - : 

C  Si  30  0  GO  TO  1400,450) ,LST 

_ a _ _ _ _ 

*TT'0 - 

576 

r77 

400 

4-tftn  ft 

WRITE  CM, 1000)  (K0UTCNP),ROUTCNP> 

FACMATI1"  f  4 

, ZOUT (NP),NP=1,NHAX) 

578 

C7< 

IwU  IT 

r*  _ 

rwn^i  ryt  trccixi  #  -  — ...  — 

GO  TO  500 

580 

Hr 

450 

- 

WRITECH, 3000)  (KGUT  CNP) ,ROUT  CNP) 
FORNATC5  (I8,SFR,  e)  ) - 

, ZOUT ( N° ) ,NP=1  »NMA  X ) 

582 

- 5 4ri - 

500  CONTINUE 


58  A 


PROGRAM  PARKTOC  /FIELD 
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SUBROUTINE  FIELD  587 


COMPUTE  ELECTRIC  FIELD  FROM  GIVEN  FIXEO  POTENTIALS  OR  CHARGES  ON  S 


1  MOMMA X, TIME »I PRINT ,R AOIUS* RRC11) ,ZZ( 11) ,X<51> • 


■j  4  J  VT  UJ  VJkftl 


tr*a  ai'i)  i 


CCMMCN/FLO/OEBYE? OIELEC  » IRELAX, PHI (11) ,0ENCC 


■Miiiik  fi  nn  i ii>  ■  l.) I 


■nil  Lirvil 


C0MM0N/CCG/N,I> J*NN, NS,NEf NH ,NSEH » CN, CS ,CE» C W,CC,C , AREA, R , Z 


REAL  KKK1*XKK2 


uo4.a>f 


0ATA  KKK1/5HALPH=/ 


.  1.14  mM  1 


OAT A  CN/5HN0RTH/ 


OATA  CE/5HEAST  / 


0IA=1 


IF(JFIRST.EQ.O)  HRITE (M, 96001  IT*  TIME 


i> 


1  1PE10«3////1X,17HC0EFFICIENT  ARRAY) 


a 


53 

l 

PROGRAM  PAfiKTOC  /FIELD 

02/li/d  0 

PAGE  14 

IF(I0H.GT.N2AI  V0LUME*AREAH< JOM) *OELZ UOM+l) 


CHARG*GCN> 


ZOLO=2 


8  00  0  FGRHATf//iX,2HZC,I2,2MI=,FlO.J/ 


2  4X,13HN0RHAL  CHARGE) 


NSEW=J0UH 


CALL  CNSEW(l) 


IF«NSEW.EQ.1J  00=ON 

-iF<NSEW«E0«2)-  0O»C)S - 

IF(OZA.GT.O)  WRITE < M,  6368)  N,I, J, 00, AREA, C 


IF(C.GT.«.)  GO  TO  150 

i  _ _ f>  _ _  , _ -  —  -  -  -  - - - 

662 

- frfrZ - 

I  ~  alph=o. 

1  _  _  _  — Tf  f  r_  rn_  ?T#  -11  _  ftp  _  ?_  FQ  -  7?  1  M  7H  A|  p  u«  A>  r|iT-?| — — - - - 

664 

- - 

X  f  *  C.  4  t  L  L  I  If  0  UK'0  tAUWlCbll’Iif  '  r  «w  r  w# 

AL«E=KKK1 

666 

- frfr? -  1 

IFCOIA.GT.fl)  WRITE (N, 7777}  AL§£, N, I , J, 00, ARE A , ALPH 

w  u  •  1 

66d 

NSEW=JOUM 


CALL  CNSEWC2J 


IFC NSEW. EO. 2)  00=0W 


PROGRAM  PARKTDC  /FIELO 


0  Z/11/8G 


PAGE  IS 


r — 

55 

|  PROGRAM  PARKTOC  /CNSEH  82/11/80  PAGE 

16  I 

■ 

1 _ _ 

SUBROUTINE  CNSEN (IGO) 

704 

1  G 

c 

'  ■  —  ■  ■  ■  ■  ■■  —  -ffL 

COEFFICIENTS  NORTH  ANO  SOUTH 

-70-* - 

706 

u 

COHHON  N,I¥,IFIRST,JFIRST,JEDGE,NR,NZA,N2,NT0T,PI, IT,IR,I2, 

-7-Q-r- - 

1  • 

2  ARE  AH  (11)  f  OELZ(ll)  9RK11)  tZCZ(ll) 

. 705 - 

710 

CUHH  UH/ FLU/ OciJTt  *  □  IfcL  EC  §  IKLLAKf  PHI  (llT'f  UCNUG 
COMHON/CCC/N,I,J,NN,NS,NE,NH,NSEH,CN,CS,CE,C«,CC,C,AREA,R,Z 

711 

** 

1  ULUHt f  tMflhb  “ 

NF C IZ# JRJ  =  JR  ♦  NRMIZ-1J 

*  1 3 

714 

'  u 

—  A 

IFtIGO.GT.l)  GO  TO  550 

rl  5 - 

7ie  . 

t  o 

AREA*  D. 

71~r - 

7ia  J 

uzjvt 

IFCZ.EQ. ff.)  GO  TO  100 

"TT  7 - V 

7  21 

IFCZ.LT.Q.)  AREA=OIELEC*AREAHIJ> 

721 

722 

c 

bO  IU  dQO 

TCJ“ 

724 

■  •  ■  "X  (iff 

Ir  1  N3Cn  tuTlJ — AKUA^WCCArtrifr  "  "  ""  - 

IF (NSEN*  EQ. 2J  AREA=OIELEC*AREAH{ J> 

rg? -  | 

726  1 

j  280 

IF(NSEN«E0.1)  GO  TO  380 

rn  rrt  »a  J.IIH  ..... 

r  c  r -  1 

728  f 

i  ,  -  ^ 

RETURN 

1, 

730  [S 

1  *00 

NN=fl 

TTI -  « 

732  1 

NN=NFfI-l,J) 

J  t  . 

734  1 

£ 

r  _ p  .  . 

Ui*it  U"ll  '  £  1 

60  TO  500 

W  is  H 

736  ‘ 

TTT 

|  400 

NS=0 

i. 

713 

1  ■*  Q 

f 

NS*NFfI»l,J> 

"P»7 - -  -  | 

740 

TL4 

n 

i"*  «  A.  /»-!  m  /AY  _  _ 

742 

“VII 

TTVT7  Vt  WFTJC  •  -  ■ 

return 

744 

c 

ee* 

AAklTTkitir  _  _ 

f  V7  - - - 

74c 

Ti.  T 

c 

uun  r  uiuc  ..  ... 

- ABfU  fl 

74  a 

ORsOt 

750 

-7C  t . 

IF  (Z*GT*  Q •)  AREA^OELZCI) 

Ten_i  t  •  t  locA^nrei  er*«ei  TitA4_v 

752 

- 

s 

/ 

754 
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PROGRAM  PARKTOC  /CNSEM 
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IF(NSEM.EQ.l)  GO  TO  600 


RETURN 


IFIJ.EQ.NR)  RETURN 


OR*RRfJ*l)  - 


_  M 


GO  TO  800 

765  | 

G 

TOO 

NH*0 

- Ttl  t  4  % - AJf  A—*  _ - - - - - - - 

T66  1 

76  7  J 

IFCJ.EQ.l)  RETURN 


0R=R  -  RR(J-l) 


IFtOR.GT.O.)  C*AREA/OR 


END 


.4  4ll|,j. 


PROGRAM  PARKTOC  /ARRAYS 


57 


02/11/8 Q 


PAGE  18 


SUBROUTINE  ARRAYS 


c 

SAVE  COEFFICIENTS  ANO  INOICES  IN  ARRAYS 

780 

w 

COMMON  MyIVylFIRST, JFIRSTy JE0GEyNR»NZAyN2*NTCT yPIylTyIRylZy 

-f*l - 

* 

1  nUIIH  f  J*r  *IMT  |RA  DlUS  1X1  yZZCXli  fXC91)  oQt/Di 

2  AREAH(ll)  yOELZ(ll)  yRIIll)  yZCI(ll) 

w- 

784 

tCniUJfJI/FLU  /  UI LL Lty  IRELIXf  PHI  f  111  | OcNGG 

COMHCN/ARRAYC/  INDXC5iy4) ,C0EF(51,5) 

r&5 - 

76  6 

- e — 

COnrON/vGu/ Hfl  >  «J  f  NWf  N  J  f  NEf  CNf  G3f  Gcf  wv*j(C#U$  j  Z 

1  y VOLUME yCHARG 

788 

-w - 

C  790 


C  810 


VVU  J  1  ^  f  —  w  s  TW  » — A  VI  I  F  r\ «/  yv  -  w 

CS8  IFIOE8YE.GT.O.)  COEFf N,5)*CHCN)/CC/0E8YE**2*.9  812 


C  CM=NET  POSITIVE  CHARGE  IN  N-TH  BOX  IF  ACT  OR  .9  FOR  CH  IN  FICOCCILOM 


IF( JFIRST.EC.O)  WRITE <M, 1000)  HN,  CN,N W, CH, N , CC , NEy CE» NS,CS, CHA ?G, 


1009  FORMAT (/13X,5I1X,1HI,I4»  2H)  =»  1PE10 «4)  ,2Ei4.4l 

- RETURN - 

END 


18 


program  parktoc  /relax 


0  2/11/8  0 


PAGE  19 


SUBROUTINE  RELAX 

T - * - ' - 

C  POINT-SUCCESSIVE  RELAXATION  TO  SOLVE  EQUATIONS 

x - - - - - - - ' - ' - - - ’ - 

COHHON  H*  IV flFIRST* JFIRS T»  JEOGE*  NR»NZA»  N2»NT0T *PI» IT* IR* 12* 

- f-RONHAX>TTWe,IPRINT>Rt>BIUS>RR<ll)  t E7t 11> » X t 5 1>»8 15 II , - 

2  AREAH(ll) yOELZfll) ,RI till ,ZCIC11) 

- COWHOW/PLO/OgSYgyOrgLtCyXReLAXyPHKll^OgNCC - 

COHMON/ARRAYC/  IN0X(51,4) ,C0EF<51,51 

X - - - - - - - : - - - 

0M€GA*1.9 

- gySjU.g,; - - - - - — - - - 

ITNAX*2fl00 

- ITRar« - - - - - - - - - - 

IPROLO=0 

- rtjt3art - - - - - 

IF { JFIRST . GT .0 .AND . IT* GT  *0 )  GO  TO  200 

- PO-Htf-  <-l,NTOT - - — - - 

100  X<K)*0. 


200  CONTINUE 

- - ITR-ITR»1 - — - - - ; - - - — 

OELTAMO* 

X - -  . ■■ — — « — = - - - - 

00  SOO  N-l,NTOT 

- X1=XIK>  - — . . g -  -----  :  - - - 

I*  ( N-ll/NR  +  1 

- JaHOOtN»t«NR)»l - - - -  — * - 

IFIJ.GT.JEOGE.OR.ZZCI) .NE.Q.)  GO  TO  NOR 

X - * - - - -f"  .  M--  I  - - - - 

C  SET  SURFACE  POTENTIALS 

-0 - — - -  . . . . — - 

X(NI*PHI<J> 

- GO-  TO  900 - - - - - - — - 

40  0  CONTINUE 

X - — - -  — - — . . — — — 

SUM*COEFlN,51 

- pp  380  K**l  *+ - - - - - 

INDEX* INOX (N,KK) 

- IFt INDEX ♦GT«0)  SUtMSUH  «  COgFCN*  lg*l  "XtlNPEXI 

300  CONTINUE 

X - 

XtNIsSUN 

- X<N)*OWg  Q4*XYN>-  »"  (  It  "OMEGA)  »X1 - 

DELTA* ADS (X  <NJ -Xll 

- IFtXlrNgrO.7  0etTR»A*3m<N7-XlT  ft  11 - 

IFCOELTA.’GT.OELTAM)  0ELTAM*0ELTA 

X - 

500  CONTINUE 

-* - - - - - - - - - - 

-4 

IFtITR.GT.ITHAX)  WRITE <H ,8883)  ITR 


829 

T5T 

831 

■VTf 

*54- 

83S 

***- 

837 

xx«- 

639 

*48“ 

841 

X4T 

843 

*44- 

845 

*4*“ 

847 

X4X 

849 

•w 

851 

-*5X- 

853 

*54 

855 

•**"«» 

058 

857 

-85-fr 

855 

-84*- 

861 

*€rE- 

863 

*44- 

865 
-86 *- 
867 
**8- 

866 
47*- 
871 
474- 
873 
-874- 
875 
47*- 
877 
47*- 
879 
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0  2/11/6  0  PAGE  23 


600  IF(OELTAM.GT.EPS)  GO  TO  200 

-e - - - 

C  ITERATION  CONVERGED.  PRINT  ANO  EXIT. 


700  I GO* 2  692 


- irurnin  i  uu — ro — long - - - - t>  T& 

BOO  NFFP=<NTOT/300)  ♦  1  854 


uw  Ar«wu'-A|nrrT-  - 

WRITE! M, 7777)  IT, TIME ,ITR, EPS, OELTAM, OMEGA  896 


wn%  b  i  ^  x  nv9-#  - - - 5  7; 

900  CONTINUE  898 

-*7-7*  FORMAT  (1H1/1 BHftFIELft-  eAlGUiATiflNr5Xy»HlT  »,  1 3, 5X  >9H*T-TIMe  - 

1  1PE10.3//15H  SOLUTION  AFTER, 16, 2 X, 25 HITERA TIONS  WITH  TOLERANCE, 

^  J*  ^  ^  ^1  ^  *  A  v«  i  Ol  n.  *  w  _ ^  A  «u  <»  fl  p  II  U 


■-  1  wn  f  AMvn  w*r  '^r\4.  nw  Uf  »  (Vf  j  o  rwrc.  j  rv  I  j  - 

3  1HN, 2X, 4HQ (N) ,7X,4HX (N)//>  902 


PROGRAM  PARKTOC  /IKTERF 


0  2/11/8  0  PAGE  21 


SUBROUTINE  INT£RPCR,Z, JG,IG,PHIC, RI,ZI,IR»IZ»I*T> _ 


INTERPOLATE  IN  R-Z  GRIO  TO  PINO  POTENTIAL  AND  GRADIENT 


COMHON/TT/HINT, X,Y,XOOT, TOOT, ZOOT, PHI, PH IR, PHIZ, OT 


IFCZ.EQ.  ZIC1»)  GO  TO  103 


M -VI  *4  n.  1-4 


IG-IZ-I+2  926 


AT  I  >  t  L  A  1  A  w  —  f  WW  I  V  AWW 

CONTINUE  928 


RETURN  WITHOUT  LOCATING  Z _ 9^0 


ACCEPT  IF  ZIIIG)  LESS  THAN  OR  EQUAL  TO  Z  LESS  THAN  ZI(IG-l). 


IFCZ.GE.  ZHIG-DI  GO  TO  102  984 


936 


IGsIG+1  __  _  988 


GO  TO  103  _ 940 


IFIZ.GE.ZICIG-I)  I  GO  TO  102  _  94g 


n —  * 

IG  DETERMINED  944 

_ _ _ _  - Abfi - 


CONTINUE  946 

_ _ _ _ . _ . _ — _ _ _ —A-Ar-7 - 


LOCATION  OF  fi  IN  ARRAY  948 

_ _ _ _ _  *  — 


ASSUME  Rill)  LESS  THAN  OR  EQU^L  TO  R  LESS  THAN  OR  EQUAL  TO  RKIR). 


952 


95 


15  J*2,IR  956 

-  _ _ _ ^  ^ 


IFCR.LT.RIIJ>>  GO  TO  153  95« 
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program  parktoc  /inter f 


CONTINUE 


0  2/ll/d 


°AGE  2  2 


RETURN  WITHOUT  LOCATING  R 


ACCEPT  IF  RKJG)  LESS  THAN  OR  EQUAL  TO  R  LESS  THAN  R1(JG*1) 


IF(R.GE.RI(JG*1) >  GO  TO  152 


IF(R.G£,RI(JG*1) I  GO  TO  152 


continue 


•OR.JG.LE.O)  GO  TO  999 


RETURN  IF  IG  ANO  JG  ANO  ALL  THAT  WERE  REQUESTED 


IF< INT «£Q«  Z)  GO  TO  3lQ 


SET  UP  FRACTIONS 


OEL2=ZIIIG-l)-ZI  (IG) 


F2= (Z-ZI (IG) l/OELZ 


P11*PHICCJG,  IG  ) 


02*(P21-Pil)/0EL  7  999 


04=(P12-P11)/OELR  1001 


CONTINUE 


INTERFOLATE  TO  FINO  POTENTIAL  ANO  COMPONENTS  OF  GRADIENT  1005 


PHIZ=02  ♦  FR* (01-02)  1007 

-PHIR^PL  -»  -PZUO3-04) - - IfiOG- 

PHIsPll  ♦  FR*(PlU-Pii)  ♦  FZ* (P211r Pli)  ♦  FR*FZ*  (P22-P21-P12+P11) 


c 

1012 

c 

c 

CANNOT  LOCATE  EITHER  R  OR  Z  OR  BOTH 

- rtrrs 

1014 

c 

999 

WRITE (MI NT, 9999)  R,Z,IG,JG 

" '  ""  .  . . 19X9  "" 

1016 

CALL  A3NOR ( 10M  ST0P22) 

STOF22 

1017 

1018 

9999  F0R«AT(///iX,32HCANH0T  LOCATE  *  OR  Z.  R»Z,IG,JG=,lP2Ei2.4,2I5> 


63 


OROGPAM  PARKTDC  /TRACK 


0  2/11/13 


PAGE  24 


SUBROUTINE  TRACK CR,  Z3  10  22 

______ _  _ — - - — 10-M - 

ADVANCE  X,Y,Z,XOOT,YOOT,ZOOT  BY  CONSTANT  ACCELERATION.  WITH  STEP-  C 


COHNON/TT/HINT,X,Y,XOOT, YOQT,ZDOT, PHI, PHIR, PHIZ, OT 


IF(R.E0.0.)  FHIX=0. 


IFtR.EC. 0.1  GO  TO  100 


PHIX=PHIR*X/R 


1026 


Y=Y  ♦  trr*{YOCT  -  .25*0T*PHIY» 


1032 


1C  34 

— - 

1036 

~ ±1-3-? - 

1038 


YOOT=YOOT  -  .5*DT*PHIY 


1042 


PROGRAM  PARKTOC  /ERF 


0  2/11/80 


PAGE  25 


FUNCTION  ERF  (X, ERROR) 


ERROR  FUNCTION 


OIMENSION  A €7) 


1  • 15201430E-3,. 27656720 E-3„ 4306 38E-4yl« 7724538/ 


RS=ABSTX1 


■  i  a  Mint  Hi  •  m  f  i 


S=AHIN1<S, 675,1 


1047 


XU  o 

1049 


1051 


1053 


1055 


1057 


SP=SP*RS 


jn—  w n  % 

IF(ABS(SP) .LT.l.E-91  GO  TO  350 


350  ERFC*1./SM**16 


RETURN 


1067 


UU 

1065 


1071 


1073 


65 


PROGRAM  PARKTBC  /EfiFINV 


0 2/11/80 


PAGE  26 


FUNCTION  ERFINVCR) 


INVERSE  ERF  RATIONAL  APPROXIMATION 


<  |  n  I  I  I. 


Ill  ■  I  ■  1  m-  m  m  m  m  i 


r  i>*  vrrn* 


DATA  01,02,03/1.432788,.  :8 9269,. 00 1*308/ 


ERFINV*-2. 

If  (R.OT.  1.  ."ORtRtLT .  8»  CALL  ABNOR  ( 10H - 

IF(R.GT.l..CR.R.LT.O.)  S  10P  33 


IFCR.LE. .09)  X»R*R0QTPI/2« 


ifif i  frin  f  —  i  n1 


P=C1.-R)  /2, 


ilrfl n-fl  fUL^MlVl1 1' 


IF(F.LE.ROUNC)  6  C  TO  150 


IF( RADIO .GT.0.)  T=SQRT (RADIO 


X0=1.+T* (Ol+T* (D2*T*03) ) 


IF ( XO.GT .0 • )  XP=XP-XN/XO 


IF(X.LT.ROUNC)  X=fl. 


150  CONTINUE 


IFCT.GT. 0.)  KAOIC*T  -  .5*AL0S«T> 


1077 

8-7-3 - 

1079 


1081 


1083 

-tea* - 

1085 


A  W  W 

1087 


1089 


1091 


A  V  J  i. 

1093 


1095 


1097 


1059 


1101 
-1162“ — 
1103 


1105 


PROGRAM  PARKTOC  /PFLOT 


02/11/80  PA6E  27 


SUBROUTINE  PFLOT  Cl)  *RI  ,ZI,NRI,NZI,  NRO,'NZO,  RAO  IUS,  WAX*  IT)  U10 


■  ii.  wm  a  m  ■ .  &  a  *  i 


1  REXPC122) , ISYMf 26) ,U (NRI,NZX) 


1112 


DATA  ISTM  /iHZ,lHA,lH8,lHC,lH0,lHE,iHF, 1HG, lHH, 1HI , 1HJ, lhK, 1HL » 


.i.ni.uTii.urM.Lj^n.ifnM.i.fMiMTMiiMnR.t* 


DATA  SNOT  /100 • , 10 • ,9« ,8 • »7» ,6» * 5 • 


MINT=6 


XNAX=-XMIN 


i  i|| 


NZCM=NZO-l 


iiiiMtiii  M  m  i — ■» .  r 


DO  10®  11=1, NZI 


M>fcf  •  w  r-f  T.r  f  Tmil 


XMIK=AHI  N1  (  XMN,U(JJ,  II)  ) 


i»k  mmm. 


C  LIST  SYM90L  CORRESPONDENCE 


WRITE (HI NT, 3001)  XWAX , XHIN,IT 


1  6Xf ^TWENTY  TWO  INTERVALS  ASSUMED  FOR  DATA. 


XMAX 


3  16//  24X,»INTERVAL*/8X,*SYKBOL*,5X,*FR0N*,12X,*T0*> 


3002  FCRMAT(5X,I3,2X,A1,1P2E15*5) 


30"3  FORMAT (1H1//) 


r** *** 


0RR=  (RI(NRI) -RI(l) )/FL0AT(NROM) 


UW  k  W  W  k  A  -  t-  J  w 

REXPf  II)  =RE  XF  <11-11  *0Rfi 


MM  Ad I f:iaiir4*i1(  «ra  -*■  j  « <  j 


Z=ZI(l)*OZZ 


00  350  KK=1,N70 


1118 


1124 


1126 


1130 


1134 


1136 


1140 


1142 


1144 


1146 


1152 


1154 


1156 


1158 
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program  parktdc  /pflot 

82/11/80 

PAGE  28 

XF<Z*LE«  ZI (I2CUR) )  IFLG*1 

1161 

J*c  UR**  r 

INTafl 

- ttfr2- 

1163 

00  30(1  11*1  f  NkO 

X=R=REXP(II) 

tie  4 

1165 

- 

IF(X.GE.RI ( JRCUR) )  JPLG*1 

1166  ~ 
1167 

LALL  X  M  t“r  tkf  C  f  J9f  1U  9  U  J  fCX'  f  £X  §  “  X  |VfZlf  I  NT  7 
INTal 

11  U  C 

1169 

11 1?" 
c 

DETERMINE  SOX  NUMBER  FROM  BNDY  MATRIX 

-  trr  fl 

1171 

DO  258  JJ>1,Z2 

- XT*  C - 

1173 

Icum^vJ 

TB0X=(  C9N0Y  (JJ)-PHD*  IPHI-9N0Y  (JJ+l))  > 

- Ttrrn — 

1179 

Ir  ( X90X* EG#  ZrJ  ILINe  ( XIi-*iSYN(2E> 

IF (T90X. 6T. C. )  GO  TO  750 

- lire — 

1177 

75  0 

CCNriNUc 

CONTINUE 

-  XT  7  “ 

1179 

aincuiz-in  - - 

IFCNOOflSOX,?)  ♦EQ.l?  XLZtECII)slSYN(N0?Ua0X-lv 26)  ♦  !> 

—  -  llctj 

1181 

- C - 

C  MRK  CORNERS 

_ A _ 

no? 

1183 

44  Ak 

■J - 

IF4IFLG.  cO.l.AND.JFLG.EQ.1)  ILINE f  II)  «1H* 

1185 

- o 

C  H4RK  DISC 

IX  0  c 

1187 

IF  (  ZI  (IZCURI  .EQ.O..A.REXPCID.LE.RAOIUS.A. 

IFLG.EQ.l) 

11  flo 

ILINE(II)=1H- 

c 

*r  (vFvGv lq v  1)  jhouw^- jmjuw •  i 

1191 

L.l_Qa _ - 

- JV  J  TVTTTITT^C 

C 

1193 

- 

C  HAVE  we  CROSSED  A  Z  80UN0ARY????? 

1199 

_ 11  qf _ 

-t*  - 

IF(IFLG.EQ.l)  GO  TO  310 

1197 

...  i-VQA  - 

- 0 - - 

C  ELSE 

_ U^y^g  1  MftlT _ VRft-Lt _ ATI  T  MCI  9\ _ _ ilAAi _ 

1199 

l  Pfl<v  - 

3304 

MJVA  IU  tliilYI  f  UkAMblll  f  ^rvVTf 

FORMAT  (9X,12|ZA1J 

1201 

i-ana - 

C 

74  ft 

1203 

_ 4-p-ft-U-  - 

41  J 

7M  AC 

WWTTTiTTwG 

MRITE<MINT,3005)  ZHIZCURI  ♦  (ILINE(I>  ,  1=1, NRO) 

14  V _ P5 _ 3 _ A  V _ ifli  A  \ _ 

12C5 

f- 

w 

"'Ww 

7C  A 

IZCUR*IZCUR+1 

1207 

_ jj.AJ. 

. 

RETURN 

-£Nfl - 

1209 

- 1*10 - 

PROGRAM  PARKTOC  /CENSTY 


SUBROUTINE  CENSTY (INIT) 


02/11/8C 
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1211 


COMPUTE  FLUXES  ANO  CHARGE  (DENSITIES  1213 

- - - trrr  - 

COMMON  M,IV,IFIRST,4FIRST, 4E0GE, NR,NZA, NZ, NTOT,PI, IT, IR, IZ, 
r  HONH«X,TlHg,IPRINT,RAPIUS,*Rtlll  ,ZZtll)  »U(S1>  ,0(91), - 1216  ■ 

2  AREAHtll) ,OELZ(ll),RI (11) ,ZCIUi>  1217 

CCWHON/OEN/RPT  ,  ZPT  ,  AL 1  ,8El,gV,  TV0LT3,  DEN 37,  H?TS,SPtgOOTCW*i - 

1  IPART,PARTCLC2J,0£LTA, SPEED, QOISK,OZMIN,OTSEC,NINJ,  1219 

"g-TSAyCtlirnOO)  t  TYCHCjrNW  AN,  NC£SC,NCA6S  ,  NGAGT,  NCM6SC 111  rFCMCI TTlTi 

3  MOUNT (51) , XKOUN (51) , MODE, WCL QUO  1221 

'COMMON /TT/tYIKT  ,X,Y,XOOT,  YOOT,ZOOTrPWT^PWTR,PMIZ,  PT - 1222 

DIMENSION  ENCESC (2) ,ENCA8S C2> *PHIC(51) ,WT(51> ,FHYR (51) ,PPYZ (51 ) , 


DATA  ENOESC/5H  ESCA,5HPES  /»END A8S/5HA ESOR  »5l-9ED 

nurmruy  (I-1I»J - - 

IF (I PRINT «GT«0)  WRITE ( M, 1900)  IT, PARTCL , SPEEOO ,CURR 

JV=  3-1 V 

'YF'C'IF IRST  «  Ng «  0  >  60  TO  98 - - - 

CALL  RANSET (NRAN) 

*PARTS=0 - - - - - - - - - 

HITE=ZZ(i)-Z2(NZ ) 

~**gffAg*T/gi,-»FR  (MR)»RR  I  MM) - — - ■ - 

AREA8=RR(NR)*HITE 

ROOTPlaSQRTCPI)  • - - - — — — - - - ~ 

NTOP=fl 

; ■-  -  - — - — * - - - 

NSIO=0 

_T€_t - - - - - - - — - — - 

JG=1 

ttIH-T=M - - - - - •— - - - — — . . : - 

Nl=l 

tMJB-0 - - - - - - - - - — - — — - 

9LK=1H 

Rg*lNft-Hf -  - - - -  - - - - —  - 1  - - 

rewind  jv 

IFCIT.NE.fl)  GO  TO  60 
00  60  1=1 » JEGGE 

MCM6£(I)^0 - L - 

NCHGIf 1) =0 

“CCttT'tttUE - 

GO  TO  90 

-COM-TIMtlg - : - 

RE  AC (It)  IN J A, ( ( VS AVE  ( I,  J) ,1=1 , 10) , 4=1 , INJAI 

■IffEOPtiyt)  96, f 6 - 

CONTINUE 

WRITEt-JV)  IMJA,  (  (V5AYE  <  I,  J)  f  I=l»  1 0 )  >J=1»I - 

NFARTS=VS AYE C9, INJAI ♦.! 


1225 

1226 
1227 
122  6 
122? 
12-32- 
1231 


1233 


1235 

-1226- 

1237 

-1226- 

1239 


1241 


1243 


1245 


1247 


1249 

-1250- 

1251 

1292 

1253 

-2254- 

1255 

2256- 

1257 

1256- 

1259 

2260 

1261 


r^-strrA'§r, 
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CONTINUE 

ir V INTT • ( 

NC  A8S=*0 


NCACTsQ 


1264 
-tt -w— 
1266 


1268 

-1269— 

1270 


IFIIPRINT.GT.O.AKO.IPART.EQ.U  WRITECH,1020)  CRRCJ)  ,  J=ltNF> 


IFITVQLTS.EC.O.)  60  TO  110 


IFUFRINT.GT.O)  WRITE  CM ,  1100 )  I, ZZCI1  ,  CPHIC  IJ)  ,J=J1,  J2> 
-CONTINUE - 


I1*NZA*1 


J1=J1>NR 

- jg»ja-»Nn - - - ~ - - 

00  150  J«1,NP 


IFCPIOLTS.EQ.a*  GO  TO  160 

-fr50 - PHIC6N»»U<N»/TiK)t:TS - - * - 

IFC IPRINT «  GT.O )  WRITE  (N» 1100)  I,ZZ  C  I)  ,  fPHIC  I  J>  ,  J=J1,  J2) 

160 — CONTINUE - — - - - » - r - > - 

IFfNOOE.NE.l)  GO  TO  180 


C  CALCULATE  PHYR  ANO  PHYZ  HATRIC5S 


IP*MINOtI*ltNZ) 


m=n*i 


JH=HAX8tJ-l*l> 


NP=NIJCI t JP) 

NN»RtJ(I~rJM) - 

PHtRlNI*CPHICCNP)-PHIClNI01/<RRCJP)-RRfJM)» 


1282 

- tw — 

1284 


1286 

- 128-? - 

1288 

- 1209 - 

1290 


1292 

- 1293 — 

1294 

- 1295-  — 

1296 

i  A  M  ^ 


1298 

—  1299 — 
1300 

- 1391 — 

1302 


1304 


130  6 


1308 

4  9A  n 


1310 

— 1311 - 

1312 
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IF<I.eCI.NZA.AN0.J*£Q.JEtl6E*i»  PHYRCNI  *  (PHICtMP>-PMICtNP)  )  / 

- r~  iRRtjpy  *raotus? - - rrrr 

c  131* 

- KHgHTjfiHy-j* - rrrr 

NP*NIJ (IPf J)  131? 

- PHTgC»l)a(FMtC-C1W>"PMTC  tNPI  1/fIHIM)  •TT  IIP)  1 - WTT 

IFCI«EQ»  NZA, AND,J«LE, JEOGE)  PHYZ  (N)*  CPHIC(NM)  -PMIC  ( N)  )  /  1319 

- r-reztrtrt-ZTtm - mt 

IFII.EQ,NZA.ANO,J,LE,JEOGE)  PHYALT(J)s(PHZC(M)-PHZC(NP))/  1321 

- rffltwmffn - - rw? 

170  CONTINUE  1323 

— ttU - CONTINUE - tW 

C  1325 

- IPCMPTS,  MEtOT"  GO  TO"  390 - rW 

C  1327 

T - SINGLE  TRAJECTORY - fNPTS*0? - ITET 

C  1329 

- NtPmwttfpr/m. - - - tw 

BET A*8£1*P 1/180.  1331 

- XRTTEtHi  im» . At r, 8Elr  ALPHA, BETA  ,  CY - 

CS1  GO  TO  400  1333 

- - — - — - rwN- 

C  17  AUG  78,  NOTE  THAT  SINGLE  TRAJECTORY  IS  INCOMPLETE  1335 

■€ - BEGINNING  'WERE, — TEMPORARILY" THIS  CAGE  IS  JtIGT  YERHINATEO. - 

C*****  1337 

- CALL  RENORTIGH - 3TOP460>  -  -  -w  - trss - r- - - 

STOP  400  1339 

C - INITIAL  CONMTtGNG  CP  TRAJECTORtES - - — - - 13N«- 

C  1341 

-338 - CONTINUE - ! — — t - - — - - - t3*e- 

1343 

-O INJECT-  NEW  PArRTICtEG  -  — -■  -■*-  -  - - M44- 

C  ASSUMED  FLUX  =  NINJ  IONS  PER  CYCLE  1345 

- - - - - - - — — —  ■ — * - - - 1 ONE 

C  1347 

- GO  -3TG  I*lYNIN J . . . .  . .  . b -  1348 

«ANl=RANFtNRAN>  1349 

- RftKfrgRANP  TNRAtO - - — - - 

RANJ^RANF ( NRANJ  1351 

- ftftH4=RANP  tttRAN) - -  -  - - - — - 1351 

RAN5*RANFCNRANI  1353 

- RANO-RANP  INRANT  . —  . .  . .  . .  1354 

NPARTS*NPARTS+1  1355 

- INJtalNJG+t - 135-t 

C  1357 

■C - SEtECT  YERTICAL  ANO  PERPENOICUt AR  COMPONENTS  OP-YEtOCTTY - tKt 

C  1359 

- YPERP-^SQRTY  AES<ALOOCRANtm - 13011 

C  1361 

- V V€  RT^ER  P I N Y  T  A  03  tl .  -1  r»R  A  N  *»  ► - tSfr* 

IF (RAN2, GT. • 5)  YVERT*-W<fERT  1363 
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_ _ _ 

C 

1364 

c 

- 1365" 

1366 

frt^ryTTCn i  » - 

CAL1*C0S  CALLS) 

- I3tr? — 

136  £ 

c 

3AU1^3IN  »  ALln)  *'* 

- 13  k* - 

1378 

G 

C 

• 

3SLcGT  TtJP^TOH  W’OH  “ItJt  ■ —  - 

1371 

1372 

»f  iRE A3aAi?t  cCAtl  I  rl  ARC rtiRciw*A  WtUBtirr 

IF(RAN3,LE*AAft£AS)  60  TO  360 

- 137-3' 

1374 

-  -  u 

C 

PARTICLE  HITS  SIOE.  CALCULATE  X,Y,Z  AND  PHI  CC0R9 . 

- -  1  o  f  v  - 

1376 

c 

LOC  AT*4HSI OE 

14  7  t 

1378 

8E1R=ASINCS9E1) 

1379 

1380 

\ 

\ - e- 

NSI0=NSI!5*1 

13  o  1 

1382 

- L383 - 

X=  Rif  CNR)  1384 


f=0.  1386 


C 

g»RAN9*f  HIT£>*c-l  4HB - : - - 

TA  7ft  11  II  -  -  -----  - _ 

-1387 - 

1388 

_  4.7  an 

1  c36’ 

CONTINUE 

ICC7 

1390 

1 

FARTICLE  HITS  T0P/90TT0M.  CALCULATE  X,Y,Z  AND  PHI  COORD. 

■"  XOT9Z  ’ 

1392 

-4-  7AT 

EE1R=2.*PI*SAN4 

/»«r  ISP4PI  --  .  ...  . 

14 

1394 

j-*  Q— _  _ 

T* 

tvtiwr ”  •  p  " 

S9E1*SINC9E1R> 

1396 

4  TO.7  _ 

V  • 

X=SGRT CR AN5) ♦RRCNR) 

1398 

1  7C  3 

Z-ZZ1NZ) 

TF#r9*  *  «  ▼  A  1  Jt-MH  v 

1400 

ii.fti 

IFf  INIT.  EQ.  1)  Z*HITE*R  AN6AZZ.CNZ) 

A  AAAT«£UAn7Tf»^ -  _  -  -  - 

1402 

--  ikfl.i _ 

L UW R  1  "UIIJU  *  1  v««  . 

IFCCALl.LT. 0.)  L0CAT*3HT0P 

1404 
■  iLfl  C _ 

IFCCALl.LT. 0.)  NT0P*NT0PA1 

FflNTTNIlF  -  -  - 

140  6 

Aun?  . 

t 

--  **  i*A 

•  AIII  ATC  an  AATT1I  AAMflAUtklte  - 

1408 

-  - 

C 

1410 
_ 1  U 1-1  - 

: 

YOOT*SP£EDO*VPERP*S9E1 

mAT.Mcen.HMCB* 

1412 

1411  ... 

C 

1414 
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1415 


VSAVE t3, IN J8)=Z 


VSAV£(5, IN JB)=YOOT 


.  Ml  t-1  |.rf«.  wt.i  r  im 


VSAVE (7, 1NJ9)=QELTA*RAN6 


Mklin4rrt f BM  \  i4fj 


IF(IFART.EQ.l)  VSA VE < 8 , IN J9) *1 . 


iij  rirvarnriB  m  i  mi 


VSAVE(9* INJ€)=FL0AT(NPARTS>+*1 


AL1=AL1R*180./PI 


K.rAR  t.i^i  k^>.i  rmg 


WRITE  TJV)  INJ9, ( (VSAVE (K  ,  J) , K*i, 1 0 ) , J*l, INJB ) 


CSS  WRITE (H, *)  ■  AT08T0C  =" ,NT0P,NSID, NBOT 


ij.rift  m  |i4frr Mi 


1  1=1,10)  ,J=i, INJB) 


IF(INIT.EQ.l)  RETURN 


IV=3-JV 


IFCIPRINT.GT.01  WRITEfM, 3410)  IT, TIME .PARTCL 


KOUNT(I) =fl 


820  CONTINUE 


IF(IPARt.EQ.l)  IIPAR=1HI 


REWIND  IV 


RE  AO  (IV)  INJA,  ((VSAVE  (I*J),  1=1,10),  J=1,INJA) 


870  CONTINUE 


NPART=VSAVE  (9,IN  J)  ♦.! 


j.  m  l  i>4  i.riikP 


IF(IPART.£').2.ANC.VSAVE(a,XNJ).GT.0.)  GO  TO  860 


1417 


1419 


1*21 


1423 


1425 


1427 


1429 


1431 


1433 


1435 


1437 


A'TW  V 

1439 


1441 


1443 


1447 


1449 


1451 


1453 


1457 


1459 


1461 


-2.,  IMPLIES  NORMAL,  ESCAPING,  OR  ABSORBEC 
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FAYTsfl. 

PATg*lH - 

ZOL  0*  VSA  V  E  (  3  f  I N  Jl 


Y*VSAVE<2,INJ> 


'  K  ] 


RsSGRT{X*X+Y«Y> 


«»ifc  m  jki  '  | 


Y00T»¥SAVE<5,INJ) 


0T=VSAWEC7,IKJ>*C2MIN 
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1466 

-1 

1463 

a  i  f  _ 

IyvT 

1470 


1472 


1474 


1476 


1478 


IF(1NTER«  EQ.fl)  CALL  I NTERP CR, Z , JG, IG, PMIC, R« , ZZ , NR , NZ , INT) 


L •  %  w  i  ■  u/ 

E=PHI*SPEED*SPEEO 


VEL*SQRT  IXOOT*  XOOT+YOOT*YOOT*ZOOT  *  ZOO  T ) 


jl  t  lArnAni  •  v  (  f  ru  w  w  v  *  r  wn*  »  w  ^  *  >y 

1  NltU(Nl) ,X, Y,ZtR#X00T, TOOT, ZDOT, FATE, VEL 


IF(WE*.EO.fl>  GO  TO  450 

PMI A*PMI *SPEE0**e - 

IFfX.EQ.ftRCNFIl  GO  TO  420 


IFtVELSQ.GT.O.)  Z00T*$IGN*S0RT«VELSQ) , ZD0T) 


GO  TO  430 


VELSQ=XOOT**2-PHIA 


1484 


1486 

1488 


1490 

-14^1- 

1492 


1494 


1496 


149! 


IFtYELSQ.LE.O.)  XOOTa-XOOT 

Aki  YY  kll  *P 


IF ( IPRlNT »  GT.3)  WRITE<M, 2180)  NPART, I0PAR,IG , JG , 
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PHIR*PHIR-»HT (N) *PHYR< N) 


C  SIFT  FCR  SPECIAL  HALF  BOX  MARKED  BY  IFDSK 


IFU.NE.NZA.CR. J.GT.JEDGE)  GO  TO  470 
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1517 


ijnj  - 


j.i  mrtr  -i  .a .  >  a  .  iiij  j.i  aaf  »  i 


GO  TO  460 


470  CONTINUE 


460  CONTINUE 


440  CONTINUE 


l.vi 


I  VV1  ^,1 1 -^.1  Jm  ■(].■  m 


ui.n.n.t-^onp.n^tinH-yt' nn'kO « ri-A? n ' v  —  n » m 


PHYa-PHl *T VOLT  S/ TVCHG 


ERFC=1.-ERF (XPHY  ,0ENA) 


F14I,111*  iiLaTi-A* kvi  i «  «  rv «k 


I*  Tyf 


»SA¥E(10,INJI=VSAVE«10,INJI***25 


465  CONTINUE 


XDOTaXCOT/ SPEED 


■ilia r.  «  mi ^ ^ ■ 


ZOOT=ZOOT/ SPEED 


CALL  TRACK CR,ZI 


XB0T=X00T* SPEED 


ZD0T=Z00T* SPEED 


iKtrifif  BA 


VSAVEI2,INJ)=Y 


M. i  jm  pr  a 


VSAVE«4,INJ)=X00T 


kfbLkiafiriviir,i&iiai|i 


VS AYE (6* INJ)*Z00  7 


i  Ik  Ml  r  IU4  ■**' 


R=SCRT(X*X*V*Y> 


IFtR.EO. 0.)  FOQTsSQRT ( X00T*XDQT*Y 00T*Y00T) 


CHECK  WHETHER  TRAJECTORY  IS  TERMINATED 


151? 


1521 


1523 


1525 


1527 


152? 


A  7JU 

1531 


1533 


1535 


1537 


1539 


IF(MONMAX.LE.l)  VSAWE  CH,INJl*l./f  l.-SQRTCPHY/fl.»PHY>7>  1541 


1543 


* 

1545 


1547 


1549 


1551 


1553 


1555 


1557 


1559 


1561 


1563 


1565 


W  W 

1567 


75 


IF{H0DE.NE.1>  60  TO  9flO 

-e - - — 

C  GENERATE  CLOUD  ANO  DISTRIBUTE  CHARGE 


00  4S8  I -1 » NTOT 


1574 

-1575 - 

1576 

~1W - 

157? 
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W*IT£  SUMMARY  FOR  CURRENTLY  ACTIVE  PARTICLE 


,13  BMH  Mil 


LlinkRJl-V! 


.Hi  ^  i.rrji1 1’ i  —  ,  ur  i.i  nxu:i:< 


1  IG,JG,N1,U(M)  ,X»Y»Z,Rf XDOT,VD0T,ZDOT»FATE, VEL 


l.ltl  1114 


1  IG,JG,N1 jUCNl)»X,Y,Z#R»XOOT»YDOT  ,  ZDOT * FATE, VEL 


SAVE  CURRENTLY  ACTIVE  PARTICLES 


IFCFAYT.LT.fl.)  GO  TO  910 


INJB*INJ9«>1 


00  915  1*1,9 


915  CONTINUE 


910  CONTINUE 


kv  i  i 


1619 


1623 


1625 


1627 


tm.  V 

1629 


1631 


1633 


1635 


1637 


IFCINJB.GT.01  WRITEIJV1  INJ9,  I  CVSAVECI  ,J)  ,1*1,10)  ,J*1,  INJ0) 


GO  TO  890 


CONTINUE 


IFCIPRINT.EQ.8)  RETURN 


WRITEIM, 3420)  PA4TCL* C8LK, 1,1=1, JEOGE) 


1  1*1, JEOGE) 


1  1=1, JEOGE) 


N1=0 


WRITE CM, 3440)  PARTCL, t 9LK, 1,1=1, NR) 


N1=N2*1 


WRITEIM, 3430)  I,  (KOUNTCJ) , J=N1,N2) 


1641 


1643 


1645 


1647 


1645 


1653 


1655 

_ 4 


1657 


1659 


1  13H  WITH  SP6E0  =,1PE15.4,7H  CM/SEC, 4X, 22H  AN0  CURRENT  0ENSITY  *, 


C  CA78  riwwwrrrun  -  ^ 

1010  FORMAT (/  5X  ,  41H0IMENSI CNLESS  POTENTIAL  ARRAY  -  ABOVE  Z*0  1667 

- -i— 5»,3gMClN  UMTS-OF  MINUS'  rVIONS/XM-ASS) - 1 frt8 - 

2  //1X,3HR  * ,14F9 » 3/  C/4X , 14F9. 3) )  1669 
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1 


it)20  FORMAT  (/5X,  41HOIHENSIONLESS  POTENTIAL  ARRAY  -  ABOVE  7=0  1670 

- t"-SX,gOHMN -OMITS  -OF  TVI0N3) - t€rH— . 

2  //1>,3HR  -,1AF9.3/  «/4X,14F9.3>>  167? 

119  8 FORMAT  it  5  H  1rTNe,I4,gX  ,  2HZ-,-FT.  3/1  /TF16.  - 167  5 - 

1200  FORMAT  (///IX,  41H0I MEN SIGNLESS  POTENTIAL  ARRAY  -  3ELQH  Z  =  0  ) 

•130 fl - POR’HATft  R-j-t*FM3INOL(?  TRAJECTORY/ IX 7 gJHANOLES  ALFMA  ANO  OCTA  ^ - 

12F15.3,11H  DEGREES  =  ,2F15.8,8H  RADIA  NS/9H  ENERGY  = ,F15. 6 ,5H VOL TS) 

-2300 FORMAT  (1 X  *  3  t  g  191 - — —  167? - 

2100  FORMATClXtI5,iH-,Al,3I3, 1P8E11. 3 , 2X, Al ,Ell. 3 >  1673 

ggfla  FORMA TUXtlSilH-iAlrre-SflPfrEtl-raT  gX^A-l^SH — RAN-»E13 *3) - 16?'3 — 

2700  FCRMAT(/1X»2A6»5X,16H ENERGY  CHANGES  *»  1P2E10.2,  1600 

- t-gfrW  tWHOtE-egOlT-f  MA»  PCR-  -STERT-  ) - tOfi  t 

3400  FORMAT (/ 5X»  33H AS SORPTION  SUMMARY  FOR  IT  CYCLE  *,I5,  1682 

- lr-g-X»-9HAT-F-IM€  ^»lP-etfrr3T2XT2A5igRf  5X»  83HTOT Afc-  MUHOCR  A9SOR9£e-^r- 

2  I7/8X,27HNUMf'ER  A9S0RBE0  THIS  STEP  =  ,  15 ,8X ,  15FNUM 8ER  ESCAPING, 

- 3-1-gH  THIS  STEP  ,  6'X  ,g5HNtfM2SR'~  CURRENTLY-  ACT  1VE  =,191- - 1665 - 

341*  FORMAT (1H1//,5X, 30 HPO PUL AT ION  SUMMARY  FOR  ITERATION  IT  =, 

- f I-S  »9X-y9MAT-TtMe  ^TtPgl* .  3f2Xf2A*>) - l~6fr* - 

3415  FORMATC5X  , 6HCELT  As , IP El 1*3  ,5X»11HC MEANS  OT  =,E11.3,5H  SEC) 

- t-2/10X-,FHl  J - Nt 2X~r9+HPOT €NTl~At."fO * y  1H X «  1 8 X>  lHYflQX  ,  1H? rt 3 ?» 1HR-, - 

2  8X,42HXOOT  YOOT  Z90T  VELOCITY,  1690 

- 3-/-g*XxgH¥0LT?»9*t:gHCn»  9X  fgHCHy  OX»  gHCH  y  9X  »  gf)CM»  7X  ,6  HCM/-SEC-, - 

4  6Xf6FCM/3EC,4X,6HCN/SEC,6X,6HCM/SEC)  1692 

■■■34  38"  FORMAT  (//  5 Xr3iHSUMHAR¥  Or  AQG0R9ED  CHANGES- FOR,2Xy  &A5//12X,9  (A-2-, - 

1  2HR(,I1,1H)),11(A1,2HR(,I2»1H)))  1694 


3440  FORMAT (//5X , 34HSUMMARY  OF  SPACE  CHARGE  MATRIX  FOR, 2X , 2 A5 // 1 2X , 9 (A  2 


1698 


p  «-■ 
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SUBROUTINE  CLOUO (RCO, ZCD,PHIC, IFDSK, MT » WCLOUO)  1699 

n++***  i7oi 

C - DISTRIBUTION  OF  FINIT5«CL0W  PARTICLE'  AHONC  GRIP  BOXES - tTGT 

C  (CHARGES  ANO  FORCES)  1703 

-c - • - - - trttr 

C  VERSION  OF  15  JAN  79.  LAST  ALTERATION  2  MAR  79.....  1705 

Z~++~* - - - — - ffOf 

C*»  ***»»•*-•****  *»*..****.**....♦***»****.***♦.***♦**..♦**.********♦** 


~C 

c 


r. 


- OIHENSION  -FWTCt  5TTl  NTtSl)  ,  IFOSKCl  t) - 17178" 

COMMON  N,IV,IFIRST,JFIRST, JEOGE,NR ,NZA , NZ, NTOT, PI, IT,IR,IZ, 

— 1  M  O  MW  A  X,  TfWiiTP  RlNT  >  R  tPTWfr  RR  ( 1 1)  f!?Ctt>Tm!irtO(>m - tTtlT 

2  AREAH(ll)  ,OELZCll)  ,RltUJ  ,ZCHil>  1711 

— nt  j  tr»  jr^HEtr  “grr.r--1  - — - — — - - trrr 

NZAP=NZA»i  1713 

- __ - - - - - tm 

DETERMINE  R  ANO  Z  CCOROINATES  OF  CORNERS  OF  CLOUO  1715 

- - - _ - tfrt 

ZH=AMINl(ZCO*WCLOUOfZCI(ll>  1717 

- ZL-  AM  A  XI  (  ZeO-W  CL'OUOyZ  Cl  MZ1 1 - t7T"8 

RH=AMIN1(RC0-»HCL0U0,RI(IR))  1719 

- RL=AHAX1  (RCC"NCLOUOr»S-) - : - 1T?8~ 

17  21 

OETgRMlMg  Ift-AKC  JO  CORNERS  OF  BRIO  CO  NT  AIMED  WITHIN  CLOW - tf«- 

1723 

- CAt fc~ TNTg-RP-fRt.  ,  2t»  JSLylSt  »~PHIC  »RI>ZCI  »TR»I2  f  23 - 1774 

CALL  INTERP CfiH,ZH,JGH, IGH,PHIC,RI »ZCI»IR*IZ» 2)  1725 

- — — ■  . . - - - - - - trftfr 


DETERMINE  VOLUME  WEIGHTS  OF  EACH  IG  JG  BOX  CONTAINED  IN  CLOUD 


-O - — ■ : - =rr - — - - - - - • - t7B«- 

DO  80  11=1, NTOT  1725 

- - - “7 - — - - - - - - - - : - 17W 

80  CONTINUE  1731 

-O - - - — - - - — * - - - ' — : - 1-^3-2- 

00  90  11=1, NR  1733 


- fP03*XID*=fl - ■ - 1784 

90  CONTINUE  1735 

0 - - - 1 T3fr 

00  10  1= IGH, IGL  1737 

- 00  20-  J=JGL,JGH - 1708 

N1=NIJ (I , J)  1735 

- Wg-A-REAMtJ) - tf^rfr 

IFCI.LE.NZA)  Wl=OELZ(I-l)  1741 

- ttM-TrgQ.NZAP)  HI -OgLHWZ A)  tOEtZ(MZAP) - 1747 


IFCI.GT.NZAP)  W1=0ELZ(II  1743 

- trtx-i-goii  rewi-  w  i-  zm-zci  t  igh) - rr** — 

IFCI.EQ.IGL)  W1=ZCICIGL-1)-ZL  1745 

- if  c-rsn.gQYi-Gt)  wi-g.»-wcL0t)0 - tf4f — 

C  174  7 

- —  IF-c-^v  fUr  J€  t>-  Wg=PI  *  <RI  ( J^L  +  l)  **2^RL,,2) - 1748 - 

IF(J.EQ.JGH)  H2=BI* CRH* *2*RI ( JGH)  **Z)  1749 


39 


PROGRAM  PARKTOC  /CLOUD 


79 


0  2/11/8  9  PAGE  40 


IP(JGH.EQ.JGL)  H2=2.*NCLOUO 

1750 

c 

IFC ZL.GT.O..CR.ZH.lT.fl..OR.J.GT.JEOGE)  GO  TO  Zl 

- 1~751 - 

1752 

IFTZCO.LT.O..ANO.I.LE.NZAP)  GO  TO  48 

- 175  7 - 

1754 

GO  TO  SO 

- tf55 - 

1756 

& 

c 

90X  IS  IN  SHADOW 

- rW —  -  - 

1758 

u 

50 

CONTINUE 

1760 

w 

Wi=0* 

hmpmb 

- e — 

Ir  Ntlrl  tucLl  *rtl  At  fCltr 

GO  TO  30 

40  CONTINUE  1766 


S25 


PROGRAM  PARKTOC  /AENOR 


0  2/11/6  0 


RAGE  41 


SUSROUTI NE  AENORCERR) 


1790 


XI  ~3X 

1792 


81 


PROGRAM  SUMRTE 


0  2/11/8  3  PAGE 


PROGRAM  SUMRTE (TAPE1, INPIT, OUTPUT, TAPE2,TAP£6*CUTPUT,TAPE5=INPLT) 


1  KHA8SI< 12) ,KHACTI<12) ,N£SG<5) ,OSET(4> 


r-m  M  u  i  rr«j 


NAMELIST  /CONTI/  OENCC  ,  A f*N,AMX  ,3MN, BMX,CMN,CMX, OMN  ,OMX 


C  OEFAULT  CONTL  ASSIGNMENTS 


OENCOIOQOQ. 


AMX=30. 


6MX=5000. 

- &HN=g50fli - — - 

CMX=3000. 


0NX=4. 


CALL  CONNEC (6) 


2002  FORM ATI IX,* ENTER  OENCC  A  AO  RANGES  OF  PRINTER  PLOT  VARIABLES*, 


REAO (5, CONTL) 


2004  FORMAT (IX, *ENTER  THIRTY  CHARACTER  IDENTIFIER*/) 


PROGRAM  SUMRTE 


02/11/80 


PAGE  2 


w 

I 


i 


C  INITIALIZE  MATRICES  FOR  ANALYSIS  OF  VARIANCES  61 

^  - - - " ~  ■  1  -  1  -  ■  ■  ■■  0  g 

00  620  1=1,12  63 

- *MA8sim»o - - - - - ev 

KHACTI  (11=0  65 

- *H*reg  cn-o — — - — - «r 

KHACTECII=0  67 

-62U - CONTINUE - frt 

650  CONTINUE  65 

- RgftmifSOOl)  IT,  TlMgyNIASS, Wiese,  NIACrrNPAeS,Nge-S&TNEACT-, - TIT 


1  RMIABS, RE MAES ,RMI ACT ,RME ACT , AYI ASS, A VE ASS, A VI ACT , AVE ACT, XMAX 


- JtTNggJLTNg  +  l - - - — - 

IF(M0O<JLIN£,50) .EQ.l)  WRITE (6 , 2 00 1 >  OENCC 
2001  -FORMAT  I1M1 , 9Xt  *TA~80t  ATEO  STATISTICAL  SUMMARY  OF  ORTA-  ■«-*-» 
♦•FUNCTION  OF  TIME.  OENCC  =  *,1P612.4I 

- trtMOOrJLTTTC-fye»-TgQ.l>  NRITg(6,2B05>  MESS - ' - 

2005  FORMAT (10X,5A1Q) 

- Ig'ttfOO t JL I'Ng , 99 >  <NHTtttTW« - 

JV  AR=MOO ( JL IKE , 1 0) ♦ 1 

- tFtEOf  m  F  - 

660  CONTINUE 

C - r - - - - - -r - 

C  ADO  TO  REMAINING  VARIANCE  SUMS 

-e - - - - - - - - - 

KHA0SHlll=KhA5Sllll)  *NI  ABS—KHA8SI ( JV AR) 

- ICHACT mi)  =  KHACTIMIT  ♦HI  AC T-KHACTI  C  J» AR) - - 

KHABSI<12)=KHABSI«12J+NIABS*NlABS-KHA9SICJVAfi>*KHABSirJVAR) 

- K*ftCTItlrg)»KMAeTIttg)»NIACT»NIACT-KHAttTI<JVAR)  »KMA  C  T 1 f  JV  ART 

KHA8SI CJVAR|=NIAgS 

- MltACYI  «JV  AR>«NlAeT - - - 

KHA8SE (11)  =  KHA0S £(111  ♦NE ABS-KHABSE ( JV ARI 
- KWACTE <1 17=KHACT£~fll)  »NE AOT-KWAOTEt  JV AR) - - - 


Tt 

73 


75 

re 

77 

Tt- 

79 

gg 

ei 

-eg 

83 

■Mr 

25 

-eg 


gg 

91 

-92- 


KHA8SE(12)=KHA8S6C12I ♦NEABS*NEABS-KHABSE<JVAR> *KHA9SECJVAR> 


KttACTg  <1  ?T-^KHACTgttg>  »NE  ACT+NgAST-KHACTgt-JVAR)  *KMA  CTE(JYAR) 


KHA8SE ( JVAP)*NEA8S  95 

- KHACT£< JgART-iMgACT - - - - — - - -  - 

730  CONTINUE  97 

t - - - 9g 

C  3EGIN  CALCULATION  OF  VARIANCES  59 


CELT=MINfl  f JLINE,10)  101 

- <HM-*gS=et0AT(KHA9SItll»)/0ELT - tgg 

AVIACT=FLOAT  (KHACTI<11) )/OELT  1C3 

- RMI  A  83=30*  T  (FLOAT  <*H A  8  3 1  ( 1  gTh/OE  L  T  -  AV  T  A  83*  A  V  I A  e  3  ) - Hhr 

RMIACT=SaRT<CLOAT( KHACTI C12) >/OELT- AVIACT*A V IACT>  10  5 

- AVEABSsFtOAT  (KHASSgdl) )  /BELT - Igt 

AVEACT=FLOAT  (KHACTE 111) )/OELT  107 

- RMCA83-SQRT (FLOAT (KMA  83g  ( 1 gTT/BSLT* AVg A83*A¥EAB3) - Igg 

RMEACT=SORTCFLCATCKHACTE<12))/DELT-AVEACT»AVEACT)  109 


-30-0g-  ■FOR*ATfr/t//g5X,-21HSUMH  ARY  FOR  THIS — STEP,  lfrX-y  19HT£N-S~T£P  RMS- ERRORS^ - 

1  17X,17H TEN-STEP  AVERAGES/ 22X, 4HIONS, 12X,9HELECTRONS, IX, 2 (3X, 2  (7X, 


83 


PROGRAM  SUNRTE 


02/11/80 


PAGE 


2  1HI,7X, iHE>)/lXt5HCYCLE,3X,4HTlNE,2t4X,i4HA8S  ESC  ACT), 

3  2<aX»3HAP3,SWt3HAM,g»,  3HACT,3X,3HACT)  ,2X,4hPMAX) - 

WRITE16,  3801)  IT,TIHE#NIA3S,NIESC,NIACT,NEA3S,N£ESC,NEACT, 

1  RHIA93,  RHEABSiftMIACT  »RHEACT,AyiAP3»A»gAP8»  AtI»€-Tr»^ACT-rX h 
MRITE(2, 3001)  IT,TIM£ ,NI ABS,NIESC, NXACT, NEA  8S,  N£ESC,NE ACT , 


30  01  FORMATCiX,I4,iP£10.3, IX, 13, 216, 3X, 13, 216 ,3X , 0P4FS. 1 , 3X,4F3. 1 , F8. 3) 

- SO  TO  650 - - — - - — - IIS - 

99  CONTINUE  12Q 

- pewinq  g - - - ftt — 

WRITEI6, 2003)  OENCC  122 


♦*  OF  TIME.  OENCC  *  *,1P£12.4) 

-WRTTgtfr,  20  05) -M€SS - - - - 

WRITEI6, 2006)  AHN, AHX , 8MN, BMX»CHN*CMX , CMN, CM  X 


♦  13X ,*SCAL£D  PR0M*,5X,F8.  1,»  T0*,F8.1/  128 

’^tox~r»o — ■» — lew . Roput m on*-?  i»xi-*oe at-eo  from-*t 9»» port »* — w^-ps-vij 

♦10X,*C  =  ELECTRON  POPULATION*, 9X ,*SCALE0  FROfc*»5X ,F6.1  ,*  TO* , 

♦  P6*1/10Xt*B — * — HirK-fHON  POTENTIAL  *rt±X , *3CAL€0 FROn*,5Xff  fr,-± , - 

♦  *  T0*,F8.1//)  132 

- . eON-TlNUg - * - : - 103— 

READ <2, 3001)  IT, TIME, NIA8S,NIESC,  NIACT, NEAPS, NEESC ,NE ACT ,RMIA3S, 


IF (EOF (2) )  40,50 

CONTINUE - - - - - 

IFLC-=0 

IF«NOOIIT»lO)t€0,0»  IFtS»l - - 

00  10  11*1,100 

L-INE<II>»1H - - - 

IF(IFLG.EG.1.AND.MOO(II,10).EQ.O)  LINE(II)*1H+ 

CONTINUE - - - - - - 

DO  20  11*1,4 

I  PV-<PV£C  <  I I)-05  6T  <  1 1  >  » /PSChlll)  *10  fl.  +  l. - 

IPV=MIN0 (HAXG(IPT,1),1001 


13? 

- 139 - 

140 

- 1-rl - 

142 

- 1-43 - 

144 

- 1-45 - 

146 


23  CONTINUE 


IF(IFLG.NE.l)  WRITE C6 « 30  03)  LINE 

3003  FORMAT  (10X  ,lH*-ylOOTl,  - 

3004  FORMAT C4X, 13, 3X,1H*, 10 0A1,1H*) 


40  CONTINUE 


